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1  Relevant  Findings 

Contract  FA7000- 10-2-0003  was  in  direct  response  to  Broad  Agency  Announcement  (BAA)  with 
the  Aeronautics  Laboratory  at  the  United  States  Air  Force  Academy.  The  announcement  in  the 
BAA  in  which  this  contract  was  directly  related  is  quoted  from  the  BAA  as. 

Current  research  strengths  include  several  complementary  thrusts .  Closed 
loop  flow  control  efforts  focus  on  aero-optic  and  energy  extraction,  with 
extensive  effort  in  the  development  of  automatic  control  algorithms  and 
techniques,  experimental  flow  control  methodologies  and  CFD  simulations. 

Work  under  this  contract  produced  new,  innovative  and  theoretical  methods  for  developing  con¬ 
trol  algorithms.  In  particular  artificial  neural  networks  coupled  with  direct  adaptive  control  was  a 
new  innovative  solution  for  achieving  successful  control  of  very  high  dimensional,  non  linear  dy¬ 
namical  systems.  This  control  technique  which  is  described  in  detail  throughout  this  manuscript 
was  successfully  applied  to  a  wide  variety  of  flows.  This  flow  control  approach  proved  to  success¬ 
fully  reduce  the  drag  on  a  circular  cylinder  by  decreasing  the  amount  of  energy  in  the  von  Karman 
street,  mitigate  optical  abberations  through  a  free  unstable  shear  layer,  as  well  as  regulate  and  ex¬ 
ploit  the  asymmetric  vortex  formulation  behind  an  axi-symmetric  bluff  body  at  high  incidence. 
Each  of  these  flow  control  applications  had  inherently  different  dynamics  including  periodic  vor¬ 
tex  shedding,  separated  free  unstable  shear  layers,  and  combinations  thereof,  which  demanded  a 
large  amount  of  robustness  from  a  control  design  perspectives.  Applications  and  demonstration 
of  successful  feedback  flow  control  where  shown  both  experimentally  and  computationally.  This 
manuscript  goes  into  great  detail  on  the  theoretical  approach  which  has  been  adopted  by  the  US- 
AFA  flow  control  group  and  then  details  the  applications  to  each  of  the  fluid  dynamics  problems. 
The  report  then  summarizes  the  business  portion  of  the  contractual  agreement. 

This  material  is  based  on  research  sponsored  by  the  US  Air  Force  Academy 
under  agreement  number  FA7000-10-2-0003.  The  U.S.  Government  is  authorized 
to  reproduce  and  distribute  reprints  for  Governmental  purposes  notwith¬ 
standing  any  copyright  notation  thereon. 

The  views  and  conclusions  contained  herein  are  those  of  the  authors  and 
should  not  be  interpreted  as  necessarily  representing  the  official  policies 
or  endorsements,  either  expressed  or  implied,  of  the  US  Air  Force  Academy 
or  the  U.S.  Government. 

[All  information  and  data  Herein  is]  approved  for  public  release, 
distribution  is  unlimited. 
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3  Proposed  Research 

3.1  Approach 

Classical  control  theory  is  limited  when  dealing  with  high-dimensional,  extremely  non-linear  sys¬ 
tems  such  as  flow  fields.  New  techniques  need  to  be  established  to  make  use  of  current  control 
theories,  while  also  allowing  for  a  reasonable  design  process  for  linear,  non-linear,  or  adaptive 
control  for  complex  flow  fields.  Low  dimensional  modeling  is  the  first  step  in  synthesizing  control 
algorithms.  Computational  fluid  dynamic  (CFD)  simulations,  which  are  numeric  approximations 
of  the  Navier  Stokes  equations  seen  in  Eq  74 

=  -Vp  +  V-T  +  f,  (1) 

are  utilized  to  produce  an  open  loop  forcing  parameter  space,  typically  varying  frequency  and 
amplitude  of  the  actuation  signal.  Numerical  reduction  schemes  such  as  principle  component 
analysis  or  proper  orthogonal  decomposition  (POD,SPOD,BPOD,DPOD)  are  then  used  to  greatly 
reduce  the  full  order  system  ( n )  to  a  truncated  mode  set  (m),  such  that  m  «  n  as  seen  by  the 
following  equation. 

m 

<I>(x,t)  =  £a/0)<jp/(x).  (2) 

i=\ 

The  resulting  mode  set  shows  the  decoupling  of  time  and  space,  correspondingly  modal  amplitudes 
(a,-(t))  and  spatial  modes  (<p,-(x)).  The  highest  energy/most  dominant  modes  are  retained  in  the 
truncation  so  that  the  reduced  data  set  accurately  represents  the  full  order  system.  The  control 
engineer  will  recognize  these  time  coefficients^,- (t))  as  the  internal  states  of  system. 

Typically,  a  Galerkin  projection  is  used  to  project  the  truncated  mode  set  onto  the  Navier  Stokes 
equations,  but  that  creates  numerous  modeling  errors,  such  as  unsatisfied  boundary  conditions, 
numerical  instabilities,  and  poor  implementation  of  actuation  term.  The  research  proposed,  uti¬ 
lizes  system  identification  techniques,  such  as  weighted  least  squares,  correlation  functions,  power 
spectral  density  with  impulse  responses,  neural  networks  (ANN-ARX),  networks  with  wavelet  ra¬ 
dial  basis  transfer  functions  (WNARX),  and  other  non-linear  methods  to  formulate  a  state-space 
system.  Typically,  control  designers  desire  these  system  identification  models  to  be  linear-time- 
invariant  (LTI)  state  space  systems.  This  allows  for  very  simple  control  design  procedures.  The 
problem  with  the  LTI  approximation  is  that  fluids  are  not  a  linear  system,  as  seen  from  the  Navier- 
Stokes  equations  which  govern  fluid  flow  (Eq.  74). 

A  linear  model  breaks  down  and  insufficiently  represents  the  flow  field  rendering  it  useless  for 
any  type  of  control  design.  New  methods  or  combination  of  methods  for  both  system  identification 
and  control  development,  from  non-linear  to  adaptive  control  techniques,  need  to  be  explored. 
This  research  proposes  to  investigate  wavelet  basis  networks  (WNARX)  to  demonstrate  the  full 
capability  of  identifying  complex  flow  response  for  a  range  of  open  loop  parameters.  The  WNARX 
will  represent  a  dynamic  model  which  can  simulate  off  design  flow  cases,  serve  as  reference  signal, 
and  ultimately  predict  closed  loop  behavior  for  control  design.  The  WNARX  model  uses  the  same 
network  architecture  as  a  neural  network;  the  only  difference  is  radial  basis  functions  are  used  as 
each  neuron’s  transfer  function.  This  is  shown  by, 

=  (3) 


<9u 

~di 


+  u  Vu 
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where  u  is  the  translation  of  the  wavelet,  .v  is  the  dilation,  and  is  referred  to  as  the  mother 
wavelet,  which  is  a  radial  basis  function  in  this  case.  WNARX  models  are  much  better  suited  for 
identifying  the  frequency  rich  dynamics  of  complex,  turbulent  flow  fields.  The  overall  WNARX 
model  is  given  by, 

N 

fit)  =  ^w/'F(^-(t-M(-))+c7’t  +  /o  (4) 

i=  1 

where  vv/  are  the  weights,  N  is  the  number  of  wavelet  functions,  cT  represent  the  linear  connec¬ 
tions,  and  /o  is  the  bias.  This  proposal  will  use  this  new  system  identification  method  to  formulate 
an  extremely  low  dimensional  model  based  from  CFD  simulations  and  POD/DPOD  decomposi¬ 
tions  to  accurately  predict  closed  loop  dynamics  of  a  given  flow  field.  This  model  is  then  used  to 
perform  feedback  simulations  to  condition  control  algorithms  which  can  then  be  applied  directly 
to  CFD  simulations  and  experiments.  Previously,  this  control  design  approach  was  applied  to  the 
simple  cylinder  wake  flow  fieldFagley  et  al.  [2009]  Siegel  et  al.  [2008],  Although  successful  for 
the  very  simple  flow,  we  desired  to  extend  the  approach  to  more  complex,  turbulent,  flows,  e.g. 
free,  unstable  shear  layer  ?.  Initial  results  are  promising.  The  WNARX  method  for  formulating  a 
model  has  proven  to  be  much  better  than  the  previous  ANN-ARX  model.  Neural  networks  have 
inherent  problems.  First,  there  is  that  no  straight  forward  method  exists  for  determining  the  number 
of  hidden  neurons,  number  of  layers,  or  parameters  of  the  regression  vector.  Training  relies  heavily 
on  trial  and  error  to  find  a  combination  of  parameters  that  yields  acceptable  results.  Second,  the 
convergence  of  these  networks  depends  heavily  upon  the  initialization  of  the  weighting  matrices. 
This  can  lead  to  drastically  different  results  when  training  a  single  network  with  a  given  set  of 
parameters  twice  because  of  the  initial  random  generation  of  weights.  Third,  a  properly  trained 
network  will  behave  as  a  black  box  in  which  little  mathematical/physical  insight  can  be  gained. 
And  fourth,  training  times  are  extremely  long  due  to  multimodal  error  surfaces  which  tend  to  trap 
the  solution  within  local  minima.  Wavenets  incorporate  POD  based  initialization  of  weighting  ma¬ 
trices  which  allow  for  much  higher  convergence  times;  the  set  radial  basis  functions  also  represent 
non  linear  limit  cycle  behavior  of  these  flows  which  reduces  the  local  minima  problem. 

This  research  proposal  continues  model  reduction  work  on  the  shear  layer  and  looks  to  vali¬ 
date  the  designed  control  laws  in  both  computational  and  experimental  studies.  The  reduced  order 
model  control  design  approach  proposed  is  extremely  powerful  and  can  be  applied  to  many  differ¬ 
ent  flow  fields.  The  control  development  method  is  not  limited  to  simplified  flow  fields,  but  well 
suited  for  highly  turbulent,  chaotic  flows.  Ultimately,  the  goal  of  this  research  proposal  is  to  refine 
the  method  of  model/controller  development  while  applying  current  knowledge  and  techniques  to 
different  flow  fields.  In  accompaniment  with  the  shear  layer  modeling  efforts,  I  intend  to  use  this 
model  reduction  method  for  multiple  flows  and  actuation  interaction,  such  as  cycloidal  propeller 
for  wave  extraction,  cycloidal  propeller  for  MAV  design,  shear  layer  with  blowing  and  suction  slot, 
shear  layer  with  plasma  actuation,  and  modeling  of  structure  fluid  interaction. 
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4  Research  Synopsis 

4.1  Overview 

Traditionally,  flow  fields  are  controlled  or  manipulated  via  passive  or  open  loop  active  control 
techniques.  Passive  strategies  make  use  of  modifying  aeronautical  properties  of  the  body  to  achieve 
a  desired  flow  state.  These  methods  add  no  external  energy  to  the  flow  and  are  time  invariant. 
Protrusions  or  intrusions  will  be  added  along  the  body’s  geometry  to  induce  some  type  of  flow 
change.  These  methods  yield  only  small  performance  gains  to  certain  flow  fields,  but  are  very 
easily  implemented  onto  an  aerodynamic  design.  A  few  examples  of  this  type  of  flow  control  are 
winglets,  fins,  or  dimples  on  a  golf  ball. 

The  other  type  of  flow  control  is  active  control  which  is  broken  down  into  two  subcategories: 
open  and  closed  loop  control.  Open-loop  strategies  make  use  of  actuators  to  force  the  flow  at 
a  given  frequency  and  amplitude  to  induce  some  sort  of  desirable  change  in  flow  state.  These 
actuators  are  typically  blowing  or  suction  slots,  plasma  actuators,  or  piezoelectric  actuators.  An 
example  of  this  is  blowing  high  frequency  pulses  along  the  leading  edge  of  an  airfoil  to  create 
small  coherent  vortex  structures  which  prevents  the  onset  of  separation  and  inherently  increases 
the  angle  of  attack  at  which  stalling  occurs. 

These  above  methods  have  been  exhaustively  studied  by  fluid  dynamicists.  The  focus  of  this 
paper  is  on  closed  loop  active  flow  control.  Instead  of  the  above  open-loop  strategies,  sensors  are 
used  to  feed  back  vital  flow  characteristics  such  as,  pressure,  velocity,  temperature,  density  mea¬ 
surements,  etc.  These  measurements  formulate  an  estimate  of  the  current  flow  state.  This  estimated 
state  allows  for  state  feedback  through  some  control  algorithm  which  prescribes  an  actuation  com¬ 
mand  to  produce  a  desired  flow  state  in  the  flow  field.  This  research  contains  multiple  tasks  to 
overcome  for  the  overall  success  of  these  ideas  such  as:  development  of  control  algorithms,  state 
estimators,  and  strategic  sensor  placement.  Mainly  because  fluid  dynamics  are  composed  of  highly 
complex,  non-linear,  stochastic  processes,  the  development  of  simple,  yet  robust  control  algorithms 
and  state  estimation  methods  becomes  extremely  difficult.  The  Formulation  of  low  dimensional 
models  which  accurately  and  robustly  represent  the  flow  over  a  given  forcing  parameter  space  is 
one  approach  to  developing  such  algorithms  and  is  explored  in  extensive  detail  in  the  following 
dissertation. 

Before  any  type  of  modeling  or  control  design  efforts  can  be  completed,  a  representation  of 
the  flow  field  shall  be  selected.  Flow  fields  are  represented  by  techniques  such  as  full  order  gov¬ 
erning  equations  (the  Navier-Stokes  equations),  experimental  setups,  numerical  approximations 
of  the  governing  equations  (computational  fluid  dynamics),  linearizing  governing  flow  equations 
(potential  flow  theory),  etc.  Each  of  these  methods  are  well  suited  for  flow  field  representation  and 
only  depend  upon  the  resources  at  hand  and  the  relevant  flow  features  desired  in  the  representation. 
Scaled  models  placed  in  wind  tunnels  accurately  create  desired  flow  fields  but  really  lack  with  ease 
of  measurements.  Measuring  techniques  such  as  hot  films,  particle  image  velocimetry,  Schlieren 
imagery,  laser  doppler  velocimetry,  etc.  do  serve  as  efficient  means  of  measuring  flow  characteri¬ 
zations,  but  these  techniques  are  local  measurements  either  temporally  or  spatially.  Measurement 
noise  is  also  a  significant  problem  for  subsequent  controller  design,  model  or  state  estimation  for¬ 
mulation  .  The  governing  equations  depicted  below  by  the  Navier  Stokes  equations  do  represent 
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100%  of  the  flow  physics,  but  really  lack  in  solvability. 


f  d  u 

P  -tT+u-Vu 


-Vp  +  V-T  +  f 


Only  a  handful  of  solutions  exist  to  these  equations.  The  solutions  that  do  exist  are  mostly  re¬ 
lated  to  extremely  over  simplified  flow  fields  which  are  often  of  little  interest  from  an  engineering 
standpoint.  The  Navier-Stokes  (NS)  equations  can  be  numerically  solved,  a  discipline  known  as 
computational  fluid  dynamics  (CFD).  Typically,  these  numeric  solutions  are  only  approximate  so¬ 
lutions  to  the  NS  equations,  but  they  do  provide  extremely  insightful  information  which  can  be 
used  for  analysis,  design,  and  simulation  purposes.  CFD  is  the  process  of  solving  coupled  partial 
differential  equations  by  the  means  of  finite  difference  methods,  finite  element  methods,  or  finite 
volume  methods.  The  spatial  domain  of  the  flow  field  is  discretized  into  small  cells  which  make  up 
a  volume  grid  or  mesh.  The  cells  are  irregular  shapes  (i.e.  rectangles,  triangles,  polygons,  etc.)  and 
resolved  on  different  spatial  resolutions.  Direct  solutions  of  the  NS  equations  (DNS)  for  complex 
flows  (large  Reynolds  number  flows)  is  not  computationally  viable  on  today’s  supercomputers.  The 
resolution  of  the  discretization  needed  is  proportional  to  the  cube  of  the  Reynolds  number  (oc  Re3) 
for  an  accurate  DNS  simulation.  For  a  turbulent  flow  this  could  be  on  the  order  of  1018  cells  for  a 
given  simulation.  For  laminar  flows  this  is  not  computationally  difficult  (as  Reynolds  numbers  are 
less  than  2100),  but  as  the  increase  of  Reynolds  number,  the  solution  becomes  unrealistic.  Many 
numeric  methods  to  reduce  computation  burden  while  resolving  the  desired  scales  exist  to  accu¬ 
rately  predict  turbulent,  compressible,  highly  complex  flow  fields  (i.e.  direct  numeric  simulations, 
large  eddy  simulations,  detached  eddy  simulation,  scale  invariance,  higher  order  turbulence  mod¬ 
els,  combinations  of  methods  with  filters  for  transition  areas)  Meneveau  and  Katz  [2000],  CFD 
allows  for  an  accurate  means  of  achieving  a  representation  of  the  desired  flow  field  and  is  typically 
the  most  widely  used  among  current  research  approaches.  Other  flow  representations  such  as  Euler 
equations,  potential  flow,  panel  methods  are  used  as  well. 

Feedback  flow  control  design  can  be  broken  into  two  main  categories.  Strategies  which  make 
use  of  low  dimensional  models  and  strategies  which  use  the  model  free  approach.  It  is  argued  that 
the  model  free  approach  has  less  chance  or  reaching  a  desired  performance,  but  may  be  simpler 
to  implement,  while  the  reduced  order  model  approach  is  more  time  consuming  formulating  the 
model  and  simplifies  control  derivation. 

The  model  free  approach  utilizes  adaptive  control  techniques  to  feedback  global  flow  variables 
in  an  experimental  setting  to  improve  flow  characteristics.  Becker  et  al.  [2006]  This  method  com¬ 
pletely  abandons  modeling  procedures  and  directly  closes  the  loop.  Control  laws  such  as  adaptive 
extremum  seeking  controllers  vary  open  loop  parameters  to  produce  desired  flow  states. King  et  al. 
[2004]Moeck  et  al.  [2007]  Typically,  these  controllers  lack  desired  performance,  but  do  prove  to 
be  useful  for  initial  experimental  investigations. 

The  second  entails  using  reduced  order  modeling  procedures  to  formulate  low  dimension  nu¬ 
merical  models  for  controller  development.  These  methods,  while  relatively  time  consuming  and 
numerically  intensive  are  able  to  produce  simple,  high  performing  control  algorithms.  Current 
techniques  do  have  inherent  problems  with  modeling  and  implementation  of  actuation,  as  dis¬ 
cussed  below.  The  focus  of  this  paper  will  be  on  the  latter  method  which  produces  smarter  control 
algorithms  based  on  reduced  order  models  (ROM’s). 

For  high  performance,  accurate  control  algorithms,  an  underlying  model  needs  to  be  formulated 
which  accurately  captures  the  desired  flow  features  in  the  analyzed  flow  is  needed.  Customarily, 


9 


Fagley 


KC  Engineering,  Inc. 


April  8,  2013 


FA7000- 1 0-2-0003 .Final 


low  dimensional  modeling  for  flow  control  purposes  is  a  three  step  process. Siegel  et  al.  [2008]  A 
flow  field  in  which  feedback  control  can  improve  a  flow  variable  or  achieve  a  certain  flow  state  is 
chosen.  The  first  step  involves  acquiring  flow  data  either  by  experimentation  or  numerical  proce¬ 
dures.  Selecting  the  correct  flow  variable  to  model  is  a  very  important  step  for  model  development. 
For  example,  streamwise  velocities  may  be  used  to  model  the  vortex  shedding  behind  a  cylinder, 
or  the  pressure  field  may  be  used  to  model  the  structures  of  a  shear  layer.  Once  the  experiments  or 
numeric  simulations  are  carried  out,  the  data  is  then  condensed  or  decomposed  by  common  tech¬ 
niques  such  as  the  Karhunen  -Loeve  procedure  more  commonly  known  as  the  proper  orthogonal 
decomposition  (POD)  method  (discussed  in  more  detail  in  Section  4.2.2. 1).  The  method  of  snap¬ 
shots  provides  a  powerful  tool  for  POD  which  allows  for  a  more  accurate  decomposition  based 
on  periodic  flow  behavior.Sirovich  [1987]  The  POD  method  allows  the  higher  frequency/lower 
energy  modes  to  be  neglected,  maintaining  the  higher  energy  modes,  more  dominant  modes.  This 
will  dramatically  reduced  the  order  of  the  data.  After  the  POD  procedure  the  data  will  be  in  the 
form, 

oo 

u(x,y,t)  =  £  a.j(t)<j>j(x,y).  (6) 

7=1 

The  third  and  final  step  is  developing  a  numerical  model  for  these  time  coefficients  (. aj(t )). 
Traditionally,  a  Galerkin  Model  is  formulated.  Here  the  truncated  mode  set  (<j>j(x,y))  is  projected 
back  onto  the  Navier-Stokes  equations  (#(•))  Noack  et  al.  [2004]Rempfer  [2000]  Rowley  et  al. 
[2004]Sirisup  et  al.  [2005]  Gerhard  et  al.  [2003].  For  convenience,  suppose  the  navier  stokes 
equations  are  written  as  the  non-linear  operator  #(•) 

u  =  u(x,t)  t  >  0  x  G  N.  (7) 

at 

The  spatial  modes  are  then  projected  onto  the  left  and  right  hand  sides  of  Equation  7,  such  that 

i  =  1,2....  (8) 


Where  the  set  of  basis  functions  must  meet  the  following  requirements:  (1)  The  basis  must  be 
complete  meaning  the  basis  must  span  the  entire  set  of  the  original  flow  field.  (2)  The  basis 
columns  must  be  linear  independent.  (3)  The  eigenfunctions  must  satisfy  the  boundary  conditions 
of  the  Navier-Stokes  equations  Rempfer  [2000],  This  projection  will  yield  a  set  of  ODE’s  which 
describe  the  evolution  of  time  coefficients, 


dai(t ) 
dt 


di(al,a2,  ...)• 


(9) 


These  ODE’s  will  be  quadratic  in  nature  due  to  the  convective  term  of  the  Navier  Stokes  equations. 
These  equations  have  proven  useful  for  analyzing  stability  of  flow  fields  and  developing  simple 
control  algorithms.  The  problem  is  this  set  of  equations,  while  mathematically  accurate,  is  numer¬ 
ically  unstable  because  of  the  lack  of  satisfied  boundary  conditions.  Also,  adding  the  actuation 
term  to  the  set  of  equations  in  9  tends  to  be  very  difficult.  Actuation  dynamics  are  a  dominant 
feature  of  feedback  control.  As  in  Noack  et  al.  [2004]  and  Rowley  et  al.  [2004],  an  actuation  term 
is  superimposed  on  Equation  10,  such  that 

=  $i(aha2,...)  +  £g.  (10) 
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Where  g  is  assumed  to  be  a  body  force  acting  at  the  node  e.  This  is  a  very  strong,  and  in  some 
flows,  false  assumption  to  make  that  fluid  dynamics  interact  linearly  with  body  forces  in  the  flow. 

An  alternative  approach  is  using  system  identification  to  produce  a  model  for  these  time  coeffi¬ 
cients.  System  ID  techniques  can  be  broken  into  two  main  categories  linear  and  non-linear.  Linear 
system  ID  has  many  widespread  uses  and  applications.  Linear  models  really  allow  for  a  very  nice 
analysis  of  the  underlying  system  by  the  ability  to  compute  stability,  observability,  controllability, 
robust  and  performance  margins.  Common  methods  include  but  are  not  limited  to  the  following: 
output  error  methods,  prediction  error  methods,  subspace  model  identification  methods,  AR/AR- 
MAX/ARX  model  identification  methods.  These  linear  models  really  fail  at  capturing  non-linear 
limit  cycle/periodic  behavior  which  is  very  common  in  unsteady  fluid  dynamics. 

Fluids  are  governed  (as  discussed  previously)  by  highly  non-linear  dynamics.  Thus,  linear 
system  ID  presents  excessive  modeling  errors.  The  alternative  is  adopt  non-linear  system  ID  algo¬ 
rithms  which  are  more  difficult  to  formulate/train,  slower  to  simulate,  harder  to  analyze,  but  in  the 
end  more  accurate  than  linear  techniques  at  modeling  fluid  behavior.  A  common  approach  is  to 
use  non-linear  Volterra  kernel  identification.Lucia  et  al.  [2003]  Balajewicz  et  al.  [2009]  For  time 
invariant,  non-linear,  continuous  time  systems  Volterra  system  ID  is  very  good  at  identifying  non 
linear  behavior.  To  model  the  response,  y(t),  of  a  system  due  to  an  arbitrary  input,  u(t),  a  second 
order  Volterra  series  is  formed  such  that, 

y(t)=h0+[  h\(t  -  T)u(T)dT  +  [  [  h2{t-T\,t-T2)u(T\)u(T2)dl\dT2  (11) 

Jo  Jo  Jo 

where  h\(t)  is  a  kernel  which  identifies  the  impulse  response  of  the  system  and  /12(f)  is  the 
quadratic  kernel  which  models  non-linearities.  Higher  order  series  can  be  expanded,  but  the 
identification  of  higher  order  kernels  increases  exponentially.Lucia  et  al.  [2003]  Techniques  ex¬ 
ist  to  prune  or  numerically  soothe  the  matrix  inversion  process  to  reduce  this  computational  bur- 
den.Griffith  and  G.R.  [1999]  The  inherent  problem  with  Volterra  series  is  that  they  are  strictly  input 
output  relationships,  that  is  y{t)  =  H  x  u(t).  The  model  doesn’t  identify  internal  dynamics  of  the 
system.  Once  the  input  becomes  zero  the  output  will  go  to  zero.  These  Volterra  series  do  have 
their  uses  in  aeroelastic  systems  and  non-linear  filter  design,  but  are  not  useful  in  the  application 
to  ROM  development  for  feedback  flow  control  Balajewicz  et  al.  [2009]. 

Another  alternative  to  non-linear  system  identification  is  using  neural  networks.  Network 
topology  ID  methods  are  a  capable  of  identifying  strong  non-linearities.  They  can  also  be  modified 
to  calculate  future  outputs  based  on  previous  simulated  outputs;  thus  having  the  ability  to  model 
dynamics  of  a  system.  They  are  not  limited  to  single  input-single  output  (SISO)  systems  but  are 
capable  of  simulating  multi-input  multi-output  (MIMO)  systems.  They  do  have  some  downsides. 
Training  is  extremely  sensitive,  and  tends  to  get  stuck  in  local  minima.  Techniques  do  exist  of  in¬ 
creasing  training  effectiveness,  but  backpropagation  algorithms  mainly  rely  on  increasing  number 
of  epochs  and  general  luck  in  achieving  global  minimumLarochelle  et  al.  [2009].  This  paper  will 
explore  the  technique  of  using  Artificial  Neural  Networks  -  Auto  Regressive  exogenous  Siegel 
et  al.  [2008]  (ANN-ARX)  non-linear  identification  methods  to  produce  an  extremely  low  dimen¬ 
sional  model  for  flow  state  simulation,  reference  signals  and  control  algorithm  development.  Also, 
the  ANN-ARX  system  ID  method  will  be  extended  to  include  wavelet  basis  functions  which  are 
very  well  suited  for  modeling  the  non-linear  periodic  behavior  of  certain  flow  fields. 
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4.2  Theoretical  Background 

A  closed-loop  flow  control  system  is  comprised  of  three  basic  components:  i)  a  sensor  or  array 
of  sensors  to  sense  and  estimate  the  current  flow  state,  ii)  a  controller  that  determines  a  control 
signal  to  achieve  the  desired  flow  state,  and  iii)  an  actuator  that  introduces  the  predetermined 
perturbation  into  the  flow.  Note  that  in  an  application,  the  system  may  have  multidimensional 
inputs  and  outputs,  i.e.  the  sensor  could  be  comprised  of  an  array  of  sensors,  and  the  actuator  could 
in  fact  be  a  combination  of  multiple  actuators  such  as  side  by  side  blowing  and  suction  slots.  In 
the  following  description  of  the  feedback  control  strategies,  the  wake  behind  a  circular  cylinder  is 
considered  as  an  illustrative  example. 

Model  Independent  Approach  Involves  the  introduction  of  sensors  in  the  wake  and  using  a  con¬ 
trol  law  (usually  linear)  to  produce  a  command  to  the  actuator  that  forces  the  flow.  The 
advantages  of  this  approach  are  that: 

•  No  model  of  the  flow  field  is  required  for  controller  design 

•  Direct  feedback  eliminates  the  need  for  a  state  estimator 

•  A  simple  control  law  may  be  implemented  in  an  experimental  setup  with  relative  ease 

For  the  circular  cylinder  wake  problem,  experimental  studies  have  shown  that  linear  pro¬ 
portional  feedback  control  based  on  single  sensor  feedback  is  able  to  delay  the  onset  of 
the  wake  instability,  rendering  the  wake  stable  at  Reynolds  numbers  about  20%  higher  than 
the  unforced  case.  However,  above  this  threshold,  single-sensor  feedback  may  suppress  the 
original  vortex  shedding  mode  but  destabilize  other  modes  Roussopoulos  [1993].  While  this 
approach  is  relatively  simple  to  implement  experimentally,  the  results  are  very  limited  for 
the  absolutely  unstable  cylinder  wake. 

Direct  Navier  Stokes  Approach  This  approach  is  more  structured  as  it  applies  conventional  and 
proven  model-based  control  strategies  such  as  optimal  control  theory  to  flow  control  prob¬ 
lems.  Abergel  and  Temam  [1990]  developed  conditions  for  optimality  for  a  few  simple 
applications.  However,  from  a  control  algorithm  point  of  view,  the  complexity  of  the  flow 
physics  results  in  a  control  problem  of  very  high  dimensionality.  Even  if  this  strategy  theo¬ 
retically  can  yield  optimal  results,  implementation  in  a  real  time  system  is  not  feasible  since 
it  would  require  the  solution  of  the  Navier-Stokes  equations  in  real  time. 

Low-Dimensional  Approach  Low-dimensional  modeling  is  a  vital  building  block  for  realizing  a 
structured  model-based  closed-loop  strategy  for  flow  control.  In  light  of  the  high  complexity 
of  the  control  problem,  a  practical  procedure  is  needed  to  reduce  the  order  by  capturing 
the  essential  physical  processes  in  a  low  dimensional  model.  A  commonly  used  method  to 
achieve  this  reduction  in  order  is  proper  orthogonal  decomposition  (POD).  This  method  is  an 
optimal  approach  in  that  it  will  capture  the  largest  amount  of  the  flow  energy  in  the  fewest 
modes  of  any  decomposition  of  the  flow.  POD  has  been  successfully  used  to  identify  the 
characteristic  features,  or  modes,  of  a  cylinder  wake  Gillies  [1998],  Siegel  et  al.  [2008], 

The  major  building  blocks  of  this  low-dimensional  approach  are  a  reduced-order  POD  model, 
a  state  estimator  and  a  controller.  The  desired  POD  model  contains  an  adequate  number  of 
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modes  to  enable  reasonable  modeling  of  the  temporal  and  spatial  characteristics  of  the  large 
scale  coherent  structures  inherent  in  the  flow  Holmes  et  al.  [1996], 

For  low-dimensional  control  schemes  to  be  implemented,  a  real-time  estimation  of  the  modes 
present  in  the  flow  is  necessary  since  it  is  not  possible  to  measure  the  whole  flow  field  di¬ 
rectly,  especially  not  in  real-time.  Flow  field  data  (e.g.  velocity),  provided  by  either  simu¬ 
lation  or  experiment,  is  fed  into  the  POD  procedure.  The  temporal  amplitudes  of  the  POD 
modes  are  determined  by  applying  a  least  squares  fit  of  the  instantaneous  data  to  the  spatial 
eigenfunctions.  Then,  the  estimation  of  the  low-dimensional  states  is  provided,  e.g.  by  us¬ 
ing  a  linear  stochastic  estimator  (LSE).  Sensor  measurements  may  take  the  form  of  velocity 
measurements  or  body-mounted  pressure  measurements.  This  process  leads  to  the  state  and 
measurement  equations,  required  for  design  of  the  control  system.  For  practical  applications 
it  is  desirable  to  reduce  the  information  required  for  estimation  to  a  minimum.  The  require¬ 
ment  for  the  estimation  scheme  then  is  to  behave  as  a  modal  filter  that  ‘combs  out’  the  higher 
modes.  The  main  aim  of  this  approach  is  to  thereby  circumvent  the  destabilizing  effects  of 
observation  spillover  Balas  [1978],  Spillover  has  been  the  cause  for  instability  in  the  control 
of  flexible  structures  and  modal  filtering  was  found  to  be  an  effective  remedy. 

To  provide  an  overview  of  the  feedback  flow  control  design  cycle  used  in  this  research  project. 
Figure  1  shows  the  main  building  blocks  in  the  process.  The  main  premise  of  this  control  design 
approach  relies  on  an  iterative  scheme  which  tweaks  current  methodologies  to  achieve  adequate 
flow  field  to  controller  design  time  along  with  desirable  closed  loop  performance.  The  develop¬ 
ment  started  with  building  a  database  of  flow  states  based  on  CFD  simulation  results.  First,  the 
natural  (i.e.  without  any  control  input)  flow  field  was  simulated.  Then,  a  number  of  simulations 
were  performed  where  periodic  blowing  and  suction  was  used  to  introduce  disturbances  at  a  given 
frequency  and  amplitude  into  the  flow.  The  results  of  all  these  simulations  were  analyzed  using 
Proper  Orthogonal  decomposition  (POD),  which  resulted  in  POD  spatial  modes  as  well  as  the  POD 
time  coefficients  for  each  time  step  of  all  simulations.  POD  modes  and  their  adjoint  amplitudes 
for  a  forcing  scenario  provide  important  flow/forcing  interaction  characteristics.  This  interaction 
is  then  modeled  through  the  time  coefficients  with  system  identification  techniques  outlined  in  fol¬ 
lowing  sections.  The  low  dimensional  model  is  then  verified  given  off  design  forcing  parameter 
cases.  Once  an  acceptable  model  is  formulated  adaptive  control  methodologies  can  be  applied.  As 
seen  in  Figure  1  at  each  design  point  a  possible  iteration  path  to  adjust  parameters  to  increase  mod¬ 
el/controller  performance  exist  along  paths  1  through  4.  This  flowchart  will  be  strictly  followed  an 
described  within  the  following  sections. 

4.2.1  Numeric  Simulation 

The  framework  of  this  control  design  approach  is  based  upon  the  use  of  open  loop  numeric  simu¬ 
lations.  The  simulations  were  performed  using  COBALT  from  Cobalt  Solutions,  LLC,  a  commer¬ 
cial  unstructured  finite-volume  code  developed  for  the  solution  of  the  compressible  Navier-Stokes 
equations.  The  basic  algorithm  is  described  in  Strang  et  al.  [1999]  and  Grismer  et  al.  [1998], 
although  substantial  improvements  have  been  made  since  then.  The  numerical  method  is  a  cell- 
centered  finite  volume  approach  applicable  to  arbitrary  cell  topologies  (e.g,  hexahedra,  prisms, 
tetrahedra).  The  spatial  operator  uses  the  exact  Riemann  Solver  of  ?,  least  squares  gradient  calcu¬ 
lations  using  QR  factorization  to  provide  second  order  accuracy  in  space,  and  TVD  flux  limiters  to 
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Control  System 


Figure  1:  Flowchart  of  the  feedback  control  development  process. 
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limit  extremes  at  cell  faces.  A  point  implicit  method  using  analytic  first-order  Jacobians  is  used  for 
advancement  of  the  discretized  system.  For  time-accurate  computations,  a  second  order  accurate 
method  with  Newton  sub-iterations  is  employed.  For  parallel  performance,  COBALT  utilizes  the 
domain  decomposition  library  ParMETIS  [Karypis  et  al.,  1997]  to  provide  optimal  load  balancing 
with  a  minimal  surface  interface  between  zones.  Communication  between  processors  is  achieved 
using  Message  Passing  Interface  (MPI),  with  parallel  efficiencies  above  95%  on  as  many  as  1024 
processors  Grismer  et  al.  [1998]. 

The  main  methods  for  calculating  complex  flows  with  a  CFD  solver  are  Direct  Numerical 
Simulation  (DNS),  Large  Eddy  Simulation  (LES),  and  Reynolds-averaged  Navier  Stokes  (RANS). 
The  RANS  approach  is  the  most  economical  since  it  is  designed  to  solve  for  the  mean  flow,  but  it 
is  subject  to  many  modeling  approximations.  Since  it  models  rather  than  resolves  the  bulk  if  not 
all  of  the  turbulent  motions,  this  would  be  an  inappropriate  choice  for  the  current  investigation. 
DNS,  on  the  other  hand,  makes  no  modeling  assumption  but  is  the  most  expensive  approach  since 
all  turbulent  motions  must  be  resolved  by  the  grid.  Since  the  smallest  scales  of  turbulence  (the 
Kolmogorov  length  scale)  decrease  rapidly  with  increasing  Reynolds  number,  this  approach  is 
limited  to  relatively  low  Reynolds  number  flows.  LES  is  less  expensive  than  DNS  since  it  models 
the  small  subgrid  scales  of  motion  and  resolves  the  rest  of  the  turbulent  motions.  However,  since 
the  large  scales  in  the  boundary  layer  are  on  the  order  of  the  boundary  layer  thickness  (which 
is  quite  thin  for  high  Reynolds  number  flows),  this  method  is  cost  prohibitive  at  high  Reynolds 
numbers  for  wall  bounded  flows. 

Detached-Eddy  Simulation  (DES)  is  a  by  now  well  known  hybrid  technique  [Spalart,  2000] 
for  prediction  of  turbulent  flows  at  high  Reynolds  numbers  [see  Spalart,  2000],  The  motivation  for 
developing  DES  was  to  combine  the  most  favorable  aspects  of  RANS  and  LES,  i.e.  the  acceptable 
predictions  using  RANS  models  of  thin,  near  wall  shear  layers  (e.g.  attached  boundary  layers)  and 
LES  for  resolution  of  time-dependent,  three-dimensional  large  eddies  (e.g.  free  shear  layers).  For 
natural  applications  of  DES,  RANS  is  applied  in  the  boundary  layer,  while  outside  the  boundary 
layer  in  the  separated  region,  LES  is  used.  An  array  of  flows  ranging  from  building  block  appli¬ 
cations  such  as  the  flow  over  a  cylinder,  sphere,  aircraft  forebody,  and  missile  base  to  complex 
geometries  including  full  aircraft  have  been  modeled  successfully  using  DES  Travin  et  al.  [1999], 
Squires  et  al.  [2001],  Constantinescu  et  al.  [2002],  Forsythe  et  al.  [2002],  Hansen  and  Forsythe 
[2003].  These  and  other  applications  illustrate  the  capability  of  DES  to  accurately  resolve  chaotic 
unsteady  features  in  the  separated  regions  along  with  a  rational  treatment  of  the  attached  bound¬ 
ary  layers.  Recent  DES  predictions  of  the  flow  around  complex  configurations  (all  using  Cobalt) 
include  the  massively  separated  flow  around  an  F-15E  at  65°  angle  of  attack  reported  by  Forsythe 
et  al.  [2004]  (this  simulation  was  the  first  eddy-resolving  simulation  of  flow  around  a  full  aircraft 
configuration),  transonic  shock-separated  flow  over  an  F/A-18E  by  Forsythe  and  Woodson  [2003], 
and  vortex  breakdown  on  an  F-18C  by  Morton  et  al.  [2003,  2004], 

These  highly  successful  simulations  demonstrate  the  capabilities  of  the  COBALT  CFD  solver. 
In  light  of  the  aero-optics  problem,  it  is  important  to  note  that  many  of  the  flows  investigated 
show  features  similar  to  the  ones  expected  to  play  a  significant  role  in  the  aero-optics  problem.  In 
particular,  the  fully  separated  flow  simulations  can  only  accurately  predict  flight  parameters  such 
as  lift  or  drag  in  a  time  dependent  manner  if  the  large  scale  motion  is  computed  correctly.  These 
same  large  scale  structures  play  a  significant  role  in  the  aero-optics  problem  as  well,  which  made 
the  COBALT  solver  a  very  well  suited  tool  for  the  numerical  aspects  of  this  research  project. 
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4.2.2  Reduced  Order  Modeling 

4.2.2.1  Numerical  Reduction  A  dynamical  system  can  be  written  as 

P(x,0  =/(x,t),  (12) 

where  x  is  defined  to  be  a  vector  over  finite  dimensional  state  space,  p  is  e.g.  the  fluid  density  at  a 
given  spatial  location  and  time,  t.  The  right  hand  side  of  the  equation  is  time  variant,  non-linear, 
and  described  by  the  Navier-Stokes  equations.  Unforced  and  forced  simulations  provide  a  flow 
data  base.  The  forcing  signal  is  limited  to  periodic  input  varying  frequency  and  amplitude.  The 
range  of  these  parameters  are  chosen  by  perturbations  to  the  natural  shedding  frequency  of  the  flow 
at  a  range  of  feasibly  implementable  actuation  limits.  This  parameter  space  provides  as  the  basis 
for  future  investigations.  The  corresponding  data  set  is  designated  by  the  matrix  £1  £  Rnxm,  where 
n  is  the  number  of  samples  in  time  and  m  is  the  number  of  spatial  grid  points. 

Data  reduction  schemes  such  as  proper  orthogonal  decomposition  (POD),  also  known  as  the 
Karhunen-Loeve  process,  have  been  used  successfully  to  reduce  the  data  to  manageable  size 
Berkooz  et  al.  [1993].  The  Method  of  Snapshots  Sirovich  [1987]  is  used  here  to  reduce  the  size  of 
the  correlation  matrix.  The  matrix  £1  £  Rnxm  can  be  decomposed  using  singular  value  decomposi¬ 
tion, 

Cl  =  UTV *,  (13) 

where  U  is  an  orthonormal  matrix  with  dimension  mxm,V*  is  also  an  orthonormal  matrix  with 
dimension  n  x  n,  E  is  a  diagonal  m  x  n  matrix  in  which  the  n  (because  typically  n  <  m)  singular 
values  are  arranged  in  decreasing  order  on  the  diagonal.  The  singular  values  of  £1  are  also  the 
eigenvalues  of  £1T£1.  Next,  define  Q  =  UY.  in  Equation  13,  which  yields 

£1  =  QV*.  (14) 

This  can  be  written  in  summation  form,  as  shown  in  Equation  15,  such  that  qi  is  the  ith  column  of 
Q;  likewise,  v,-  is  the  ith  column  of  V, 

m 

n  =  l >v*.  (15) 

1=1 

Equation  15  is  still  an  identity,  i.e.  no  approximations  have  been  introduced.  In  Equation  15  the 
ith  temporal  coefficient,  cij(t),  is  exactly  equivalent  to  the  ith  column  of  Q.  Likewise,  the  ith  spatial 
mode,  (pi(x,y ),  is  represented  by  the  ith  row  vector  of  V*.  The  system  £1  can  then  be  written  as, 

m 

£l  =  Y,<t)(pi(x).  (16) 

i—  1 


For  w!  <  to  the  decomposition  becomes, 

m' 

£a,-(f)9*(x)-  (17) 

(=1 

where  now  the  right  hand  side  is  an  approximation  of  £1.  Because  the  singular  values  can  be 
related  to  the  energy  of  the  modes  and  are  arranged  from  largest  to  smallest  (<7i  >  <72  >  •  •  •  > 
on),  the  dominant  spatial  and  temporal  modes  appear  first  in  the  matrices  U  and  V,  respectively. 
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Therefore,  plots  of  the  singular  values  are  normally  shown  to  determine  where  a  truncation  of  a 
model  is  appropriate.  For  example,  for  the  cylinder  wake  flow  field  98%  of  the  energy  in  the  flow 
is  contained  within  the  first  three  modes.  Reducing  m'  to  three  represents  a  significant  reduction 
in  the  model  order  even  though  the  model  still  maintains  the  dominant  flow  physics  Siegel  et  al. 
[2008], 

The  above  decomposition  is  a  good  approximation  for  a  periodic  flow  field,  but  in  feedback 
control,  the  flow  is  being  forced  through  some  (arbitrary)  actuation  and  therefore  not  periodic, 
even  if  it  would  be  naturally.  This  means  that  the  flow  will  undergo  some  transient  development, 
which  results  in  the  spatial  modes  (<p,  (x))  evolving  over  time.  Siegel  et  al.  Siegel  et  al.  [2008]  de¬ 
vised  a  strategy  called  Double  Proper  Orthogonal  Decomposition  (DPOD)  to  accurately  model  this 
transient  behavior.  In  DPOD,  the  POD  process  is  used  twice  to  capture  the  transient  phenomenon 
in  the  spatial  modes.  The  second  decomposition  represents  the  spatial  mode  fluctuations  over  time 
which  capture  forcing-flow  interaction.  The  DPOD  decomposition  is  written  as, 

m'  rl 

«= EE  «</(')«/(*)•  os) 

!=  1  ./=  1 

The  resulting  spatial  modes  form  a  m'  xn'  mode  set  which  accurately  represents  the  unforced, 
forced  and  flow  state  transitions  from  one  to  the  other.  For  more  information  on  the  DPOD  process 
seeSiegel  et  al.  [2008], 

4.2.3  System  Identification 

The  system  identification  step  as  presented  in  Fagley  et  al.  [2010]  is  a  crucial  step  in  formulating 
the  low  dimensional  model.  Numeric  decomposition  parameters  such  as:  forcing  parameters, 
spatial  domain,  and  numerical  decomposition  method  (DPOD,  SPOD,  POD.. etc.)  are  determined 
by  a  parameter  study.  The  next  task  for  defining  a  reduced  order  model  was  the  development  of  a 
mathematical  model  which  accurately  represents  the  time  coefficients  of  the  POD  mode  set.  The 
formulation  of  a  reduced  order  model  which  accurately  relates  the  forcing  input  to  the  evolution  of 
the  time  coefficients  would  give  capability  to  predict  a  flow  state  not  present  in  the  original  dataset. 
These  predicted  time  coefficients,  multiplied  with  the  spatial  modes,  would  yield  a  prediction  of 
the  flow  field  at  any  instant  in  time  within  some  confidence  interval. 

For  model  development,  the  Galerkin  method  has  been  typically  used  to  project  a  truncated 
mode  set  onto  the  Navier  Stokes  equations,  resulting  in  a  set  of  ordinary  differential  equations 
[see  e.g.  Berkooz  et  al.,  1993],  However,  it  has  been  realized  that  the  resulting  equations  are 
mathematically  unstable,  that  the  resulting  model  cannot  satisfy  the  boundary  conditions  due  to 
POD  truncation.  Also  a  linear  flow  interaction  with  actuation  is  taken  into  account  through  the 
body  force  term  in  the  Navier-Stokes  equations  Rempfer  [2000].  Nonetheless,  these  reduced  order 
models  provide  insight  into  the  flow  physics  and  basic  controller  design.  Galerkin  models  lack 
of  actuation  characterization  poorly  models  the  dynamics  of  typical  actuation  methods  such  as 
blowing  and  suction  slots  or  plasma  actuators.  The  body  force  assumption  does  not  capture  forcing 
dynamics  of  a  blowing  and  suction  slot  in  which  the  mass  flow  rate  is  actually  changing  through  a 
prescribed  location. 

These  shortcomings  of  the  Galerkin  model  drive  this  research  project  a  different  direction. 
System  identification  techniques  were  used  to  identify  a  model  which  accurately  represents  the 
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open  loop  forcing  dynamics  of  the  flow  field.  Writing  the  system  in  Equation  17  with  the  density 
of  the  flow  as  the  kernel, 


p(x,t)=f(x,t,u )  ~  £a('0,w)<jp/(x|w),  (19) 

(=i 

where  x  e  W1,  ta  <  t  <t\  and  u(t\F.A)  =  Asin(2nFt),  provides  the  basis  for  applying  system 
identification  tools.  In  order  to  understand  the  evolution  of  the  density  field,  p(x,t),  in  time,  its 
derivative  must  be  found, 


P (x, t )  ~  Y,  [di(t,u)(Pi(x\u)+ai(t,u)<Pi(x\u)]. 
i=  1 


(20) 


The  second  term,  a,-(f,  u)<pj(x\u),  drops  out  because  (p  is  time  invariant,  so  all  of  the  nonlinearities  of 
the  system  are  contained  within  the  evolution  of  the  mode  amplitudes  d,  (r|w).  System  identification 
was  then  used  to  provide  a  mathematical  representation  of  the  evolution  of  these  mode  amplitudes. 
Here  the  system  can  be  represented  in  state  space  form  as 

f  d(t,u)  =  G(a(t,u )) 

1  p(x,t)  —  (p{x\u)a{t,u) 


where  G(a(t,u ))  is  an  unknown,  non-linear,  time  invariant  function.  In  discrete  time,  this  system 
is 

f  a(tk+huk+l)  =  G(a(tk,uk)) 

1  p(x,tk)~q>(x\uk)a(tk,uk). 

The  modeling  goal  was  to  formulate  a  mathematical  model  which  represented  the  time  coefficients 
of  the  numerical  approximation  over  the  open  loop  forcing  parameter  space,  u(t\F,A).  Nonlin¬ 
ear  auto  regressive  exogenous  (NLARX)  systems  were  used  to  identify  the  behavior  of  the  mode 
amplitudes,  for  which  a  regression  vector  was  formed  such  that 


@(t)  =  [u(t) . . .  u{t  -  nu) ,  a\  (t  -  1) . . .  a\  (t  —  nai ) . . .  at(t  -  1) , . . . a,(t  -  nai)\ .  (23) 


Nonlinear  mathematical  models  were  then  trained  to  minimize  the  error  between  the  training  set 
and  predicted  output  in  a  least  squares  sense. 

Many  methods  for  nonlinear  system  identification  exist.  Previous  work  used  artificial  neural 
network  autoregressive  exogenous  (ANN-ARX)  systems  to  identify  the  dynamical  behavior  of  the 
time  coefficients  in  the  forced  cylinder  wake  Siegel  et  al.  [2008].  Neural  networks  are  widely 
used  in  the  scientific  community  for  process  modeling,  artificial  intelligence,  pattern  recognition, 
machine  learning,  etc.  This  nonlinear  system  identification  technique  has  been  argued  to  be  a 
universal  approximator,  capable  of  representing  any  type  of  data  trend  Norgaard  et  al.  [2003], 
However,  some  inherent  problems  of  ANN  models  exist.  First,  there  is  no  straight  forward  method 
for  designing  the  network,  including  determining  the  number  of  hidden  neurons,  the  number  of 
layers,  or  the  parameters  of  the  regression  vector.  Furthermore,  training  relies  heavily  on  trial  and 
error  to  find  a  combination  of  parameters  that  yields  acceptable  results.  Second,  the  convergence 
of  these  networks  depends  strongly  upon  the  initial  (usually  random)  weights  in  the  weighting 
matrices.  This  can  lead  to  drastically  different  results  when  training  a  single  network  with  different 
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sets  of  parameters.  Third,  a  properly  trained  network  will  behave  as  a  black  box  from  which  little 
mathematical  and  physical  insight  can  be  gained.  And  fourth,  training  times  are  extremely  long  due 
to  multimodal  error  surfaces  which  tend  to  trap  the  solution  in  local  minima,  which  also  contributes 
to  the  vastly  different  network  parameters  obtained  from  the  random  data  initialization. 

A  different  way  to  formulate  a  model  is  by  using  wavelet  transformations,  which  are  known 
for  their  ability  to  compress,  decompose,  and  approximate  scientific  data  sets  accurately  and  effi¬ 
ciently.  They  are  used  in  many  technical  fields  including  image  processing,  edge  detection,  large 
scale  monitoring  processes,  transient  detection,  etc.  Mathematically,  the  mother  wavelet,  VF,  can 
be  written  as 

(24) 

where  u  denotes  the  shift  or  translation  and  5  denotes  the  dilation  or  frequency  content  of  the 
wavelet  basis  function.  In  the  current  modeling  approach,  wavelets  were  used  as  transfer  functions 
to  create  a  wavelet  neural  network  (WNN).  These  wavenets  were  first  introduced  by  Zhang  et 
al.  and  have  been  applied  to  many  areas  such  as  functional  approximation,  system  identification, 
adaptive  control,  and  non-linear  modeling  and  optimization  Zhang  and  Benveniste  [1992],  Zhang 
et  al.  [1995],  Chen  and  Bruns  [1995],  M.M.  Polycarpou  and  Weaver  [1997].  An  example  of  the 
formulation  of  a  wavelet  based  ARX  network  is  Zhang  and  Benveniste  [1992] 

N 

fit )  =  Y,wiX¥isiit  -«/))+  cTt  +  fo,  (25) 

(=1 

where  Wj  are  the  weights,  VP  is  the  wavelet  function,  N  is  the  number  of  wavelet  functions,  cT 
represents  the  linear  connections,  and  /o  is  the  bias.  The  WNN  is  typically  initialized  using  a 
dyadic  wavelet  decomposition  Oussar  and  Dreyfus  [2000].  The  above  parameters  are  updated  via 
a  gradient  descent  method  to  minimize  the  cost  function 

Je  =  \\f(t)-f(t)\\2.  (26) 


Multiple  techniques  exist  to  design  the  architecture  of  such  wavenets.  One  technique  is  to  replace 
the  existing  transfer  function  of  a  neural  network  (usually  sigmoid  or  signum  functions)  with  a 
radial  basis  wavelet.  Another  approach  for  integrating  these  two  ideas  is  to  use  the  wavenet  as 
a  preprocessing  filter  for  the  non-linear  ANN  identifier.  An  example  of  this  type  of  network  is 
shown  in  Figure  2,  which  was  used  by  Angrisani  et  al.  [1998]  to  identify  transients  in  power 
signals.  This  approach  was  taken  in  the  current  research  to  design  wavenets  for  feedback  flow 
control  applications. 

The  fundamental  WN  structure  used  for  this  approach  to  model  the  system  in  Equation  22  is 


a 


n  n 

j(tk\®j)  =  (®j-r)PL+  Y,  aSif(bSi(Qj  -  r)Q-  cSl)  +  Y  awtg  (bWi(®j  ~  r)Q  -  cWi)  +d  (27) 

linear  ^  ^  s 


scaling  block 


wavelet  block 


where  r  is  the  mean  of  the  regression  vector,  P  is  the  linear  subspace,  L  are  the  linear  weights,  Q  is 
the  nonlinear  subspace,  aSi  are  the  scaling  block  coefficients,  bSi  are  the  scaling  block  dilations,  cSi 
are  the  scaling  block  translations,  aWi  are  the  wavelet  block  coefficients,  bWi  are  the  wavelet  block 
dilations,  cWj  are  the  wavelet  block  translations. 
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Figure  2:  Training  method  for  Wavenet  ARX  system  which  updates  frequency  and  ANN  parame¬ 
ters  during  single  back  propagation  iteration  Angrisani  et  al.  [1998]. 

Moreover,  f(x)  is  the  scaling  function.  Here,  the  scaling  function  was  chosen  to  be  a  class  of 
radial  basis  function  such  that 

f(x)  =e_2lWl2  3f:Rn^R.  (28) 

Likewise,  g(x)  is  the  wavelet  basis,  which  is  also  a  radial  basis  function  in  the  form 

g(x)  =  (Nlo -  N© e~^  3g:Rn^R.  (29) 

Examples  of  these  basis  functions  are  shown  in  Figure  3.  Notice  that  these  are  continuous  functions 
with  defined  derivatives  over  an  infinite  domain. 

The  linear  and  nonlinear  subspace  matrices  (P  and  Q,  respectively,  in  Equation  27)  were  initial¬ 
ized  by  a  principal  component  analysis  based  on  an  optimal  representation  of  the  system  linearities 
in  the  linear  block  as  well  as  the  nonlinear  block.  Given  a  set  of  initial  parameters  for  the  WN, 
the  model  simulation  was  performed  and  the  global  error  of  the  training  data  was  determined  as 
aj{t )  —  The  parameters  of  the  WN  were  then  updated  via  a  gradient  descent  method  to 

minimize  the  error. 

A  graphical  representation  of  this  network  is  shown  in  Figure  4.  The  regression  vector  is 
presented  to  the  three  blocks  as  discussed  above.  The  network  simulates  the  estimated  output  for 
the  entire  training  set,  computes  the  error  and  updates  the  network  parameters  in  Equation  27. 
Wavelets  as  a  set  of  basis  functions,  represented  in  Equation  27,  allow  for  a  basis  which  represents 
a  variable  frequency  domain  (by  the  adjustment  of  the  dilation  parameter,  bWi ).  The  frequency  rich, 
nonlinear  limit  cycle  behavior  of  the  two  dimensional  shear  layer  was  accurately  represented  by 
the  set  of  wavelet  basis,  [(g)(x)(/)(x)],  as  shown  below. 

A  number  of  parameters  needed  to  be  determined  to  find  a  suitable  wavenet  model  which 
accurately  represented  the  density  states  of  the  flow  field.  The  first  parameter  was  the  composition 
of  the  regression  vector,  @y.  This  regression  vector  was  composed  of  previous  time  histories  of 
simulated  mode  amplitudes  and  current  and  previous  inputs  to  the  system.  The  total  time  history 
encapsulated  by  the  regression  vector  determined  the  amount  of  memory  the  model  had.  However, 
large  regression  vectors  drastically  increase  training  times  and  may  possibly  increase  the  final 
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Figure  3:  Example  of  wavelet  transfer  functions  fix)  for  given  scaling,  dilation,  and  translation 
parameters. 


[Xj^, ...,Xn] 


Figure  4:  Graphical  representation  of  the  current  wavenet  approach  for  identifying  the  evolution 
of  time  coefficients.  Back  propagation  was  used  to  update  parameters  in  Equation  27  in  a  least 
squares  fashion. 
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simulation  error.  For  periodic  systems,  typical  time  histories  should  approximately  include  one 
fourth  of  a  period  of  the  predicted  output. 

4.2.4  Feedback  Control 

4.2.4.1  State  Estimation  The  estimator  design  process  is  extremely  important.  The  main  idea 
is  described  in  ?.  The  state  estimator  will  relate  an  array  of  surface  mounted  sensor  signals,  de¬ 
fined  as  p(xs,t),  to  the  flow  state  which  is  modeled  by  the  time  coefficients  of  a  POD  truncation 
(Oj(t)  in  equation  17)  (Note:  the  f  superscript  designates  that  the  parameter  was  inflow  and  the  ,v 
superscript  designates  that  the  parameter  was  on  surface).  The  goal  here  is  to  take  an  experimen¬ 
tally  feasible  number  of  surface  mounted  sensors  (pressure  transducers  for  example)  and  through  a 
mathematical  modeling  technique,  formulate  a  mapping  of  sensor  signals  to  the  flow  state.  Having 
access  to  the  current  flow  state  allows  for  state  feedback  flowcontrol.  The  process  to  developing 
this  mathematical  relation  is  describe  below. 


4.2.4.2  Sensor  Placement  A  heuristic  approach  to  sensor  placement  is  used  in  this  study.  Lo¬ 
cations  correlated  spatially  to  desired  flow  features  (e.g.  vortex  shedding,  vortex  pairing,  boundary 
layer  growth,  separation  points,  etc.)  are  chosen  and  defined  as  (yv)  within  the  numeric  simulation. 
A  surface  POD  analysis, 

k 

P(xs,t)  ~  J^asp(t)(Pp(xs),  (30) 

p= i 

yields  surface  POD  modes  (j>p(xs).  The  resulting  locations  of  the  maxima  and  minima  of  the  sur¬ 
face  modes  show  where  the  largest  variability  of  the  signal  occurs;  hence,  they  indicate  preferred 
locations  for  sensorsCohen  et  al.  [2003b].  The  corresponding  surface  POD  analysis  allows  for  the 
reduction  of  the  number  of  sensors  needed  to  accurately  estimate  the  surface  POD  mode  ampli¬ 
tudes. 

The  surface  time  coefficients  (a  linear  pre-filter)  are  then  computed  by  solving  for  ap(t)  in 
equation  (30),  given  a  particular  simulation,  using 

ap(t)  =  p(xs,t)(psp~l(xs).  (31) 

The  matrix  <pp(xs)  provides  the  linear  subspace,  with  dim(asp(xs ))  <  dim(xs),  on  which  the  sensor 
signals  are  projected.  The  state  vector  is  then  given  by 


p(xs,t)(pp  \xs) 
p(xs,t) 


(32) 


which  will  be  estimated  using  state  estimation  methods.  The  estimator  will  yield  a  model  for  the 
POD  time  coefficients  in  which  the  flow  state  is  estimated  by  the  linear  or  non-linear  mapping  of 
the  state  vector  through  the  function  G, 

afj(t)~G(Q(xs\t)).  (33) 

In  the  following,  the  most  prevalent  estimation  techniques  for  feedback  flow  control  are  outlined. 
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4.2.4.3  Linear  Stochastic  Estimation  Linear  stochastic  estimation  (LSE)  is  chosen  as  the 
baseline  estimation  method.  Since  flow  fields  of  interest  are  typically  highly  non-linear,  the  perfor¬ 
mance  of  this  method  usually  tends  to  be  poor.  The  LSE  method  is  also  a  static  estimator,  meaning 
no  time  histories  of  the  sensor  signals  are  used  to  improve  the  mapping  performance.  Nevertheless, 
a  linear  analysis  is  always  important  because  it  serves  as  a  benchmark  comparison  for  the  more 
complex,  non-linear  system  ID  methods.  For  the  LSE  method,  the  estimated  state  aj(t)  is  obtained 
from  a  linear  mapping  through  matrix  L  where  L  £  Rnxm  with  m  is  the  dimension  of  the  sensor 
measurements,  n  the  dimension  of  the  state  space.  The  LSE  operator  is  given  by, 

Oj(t)  =Ld(xs\t).  (34) 

The  observation  matrix,  L,  is  obtained  by  correlation  of  the  data,  such  that 


LijE(0(xs\t)i0(xs\t)j)  =E(a{  (t)d(xs\t)j), 


(35) 


where  E(.)  is  the  expected  value.  This  can  also  be  written  as, 

L  =  BA~X  where  B  =  E(a^(t)d(xs\t)T)  and  A  =  E(6(xs\t)6(xs\t)T)1  (36) 

Noting  that  the  span  of  L  is  limited  to  the  number  of  sensors.  Likewise,  the  number  of  states  must 
be  smaller  than  the  number  of  sensors  for  this  estimation  method.  This  method  allows  for  very 
quick  computation  and  simple  set  up  within  simulation  and  experimental  settings.  LSE  is  highly 
sensitive  to  noise  and  requires  a  large  array  of  sensors  to  provide  accurate  estimates. 


4.2.4.4  Artificial  Neural  Network  Estimation  Artificial  neural  network  estimators  (ANNE) 
are  a  powerful  non-linear  system  identification  methodCohen  et  al.  [2007],  A  two  layer  feed  for¬ 
ward  network  is  used  in  this  study  to  map  the  sensor  signals  to  the  current  flow  state.  The  first  layer 
implements  the  tank  transfer  function  while  the  second  layer  consists  of  linear  neurons.  Previous 
time  histories  of  the  sensor  data  are  used  in  the  ANNE  effort  such  that  the  estimation  is  autore¬ 
gressive  (AR).  The  regression  vector  x  is  formulated  by  concatenating  current  and  previous  state 
vectors, 


*0) 


e(xs,t) 

d(xs,t-ti) 


(37) 


0(xs,t-tn) 


The  regression  vector  is  then  presented  to  the  network  which  can  be  written  as 


Oj(t)  =  WotanhffixW+bti+bo, 


(38) 


where  Wj,  Wo  are  the  input  and  output  matrices,  respectively,  and  bi  and  bo  are  the  input  and  output 
biases.  During  network  training,  the  global  error  for  a  given  training  data  set  is  estimated  and  the 
weighting  matrices  along  with  the  bias  weights  are  updated  via  the  gradient  descent  method  to 
minimize  the  estimation  error. 
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4.2A.5  Wavenet  Estimation  Wavenet  estimation  (WNE)  methods  combine  the  network  archi¬ 
tecture  of  the  above  ANNE  method  with  wavelet  basis  functions.  Wavenets  were  first  introduced 
by  Zhang  et  al.  and  have  been  applied  to  many  areas  such  as  functional  approximation,  system 
identification,  adaptive  control,  and  non-linear  modeling  and  optimization.Zhang  and  Benveniste 
[1992],  Zhang  et  al.  [1995],  Chen  and  Bruns  [1995],  M.M.  Polycarpou  and  Weaver  [1997]  Multiple 
techniques  exist  to  design  the  architecture  of  such  wavenets.  One  technique  is  to  replace  the  trans¬ 
fer  functions  of  a  neural  network  with  a  wavelet  basis  function.  Another  approach  for  integrating 
these  two  ideas  is  to  use  the  wavenet  as  a  preprocessing  filter  for  the  non-linear  artificial  neural  net 
(ANN)  identifier.  In  this  research,  a  combination  of  these  two  methods  is  used.  The  model  struc¬ 
ture  is  decomposed  into  three  blocks,  a  linear  block,  a  preprocessing  scaling  block,  and  a  wavelet 
block.  The  AR  model  is  then  trained  to  accurately  estimate  the  frequency  rich,  highly  non-linear 
POD  modal  amplitudes.  The  fundamental  WN  structure  to  model  the  system  aj(t)  =  G(x(t))  is 

n  n 

afj(t)  =  (x-r)PL+Y,  aSif(bSi ( x-r)Q  -  cSi )  +  £  awtg (bWi (x-  r)Q-  cWi)  +d,  (39) 

i=  1  i=l 

where  r  is  the  mean  regression  vector,  P  is  the  linear  subspace,  Q  is  the  non  linear  subspace,  L 
are  the  linear  coefficients,  ani  are  the  wavelet  coefficients,  b,v  are  the  wavelet  dilations,  c,v  are  the 
wavelet  translations,  and  f(x)  is  the  scaling  function.  In  this  investigation,  the  scaling  functions 
were  chosen  to  be  a  class  of  radial  basis  function  such  that 

f(x)  =  g4lWl2  3f:Rn^R.  (40) 

Likewise,  g(x),  the  wavelet  basis,  is  also  a  radial  basis  function  of  the  form 

g(x)  =  (||x||0  -  \\x\\22)  e-'M  3  g  :  IT  E.  (41) 

The  linear  and  non-linear  subspace  matrices  (P  and  Q,  respectively,  equation  39)  are  initialized 
by  a  principal  component  analysis  based  on  an  optimal  representation  of  the  system  linearities 
in  the  linear  block  as  well  as  the  non  linearity  block.  Given  a  set  of  initial  parameters  for  the 
WN,  the  model  simulation  is  performed  and  the  global  error  of  the  training  data  is  determined 
as  aj(t)  —  dj(t).  The  parameters  of  the  WN  are  then  updated  via  a  gradient  descent  method  to 
minimize  the  cost  function 

Je  =  \\afj(t)-aj(t)\\2.  (42) 

4.2.4.6  Adaptive  Control  A  feedback  law  needs  to  be  developed  to  control  the  model  as  seen 
in  (99).  This  equation  is  written  again  as  follows, 

f  d(t,u)  =  G(a(t,u)) 

[  p(x,t)  ~  (p(x\u)a(t,u) 

The  WNARX  model  estimates  a  given  flow  state  based  on  forcing  input  and  dynamics  of  the  system 
as  seen  by  (27).  The  closed  loop  goal  is  to  regulate  (i.e.  force  a  particular  state  to  zero);  thus, 
reducing  active  vortex  shedding  and  achieving  a  desired  closed  loop  flow  state.  Direct  adaptive 
control  law  was  chosen  to  close  the  loop.  This  method  of  control  allows  for  unforseen  uncertainties 
when  scaling  the  controller  up  to  feedback  CFD  simulations  and  experiments. 
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The  ROM  in  (27)  allows  for  relatively  short  simulation  times  to  pre-condition  an  adaptive 
controller.  Adaptive  control  theory  demands  a  linear  system  to  prove  stability.  The  model  in  (27)is 
highly  non-linear.  Linearizing  would  lose  many  of  the  flow  effects  which  were  desired  in  the  model 
in  the  first  place.  To  develop  an  adaptive  strategy,  satisfy  stability  and  bounded  trajectory  issues, 
the  model  can  be  written  as, 

{ <«> 
Even  though  flow  fields  are  highly  non-linear  and  linearization  is  often  a  poor  approximation,  there 
is  still  information  to  be  learned  from  the  linear  model.  Suppose  the  adaptive  control  goal  is  to  have 
the  output  of  the  plant  go  to  zero,  commonly  known  as  output  regulation,  that  is 


7  >  °o 


0. 


(45) 


Suppose  there  exists  a  G*  which  moves  the  system  along  some  ideal  trajectory,  such  that  the  closed 
loop  system  can  be  written  as 

x=(A+5G*)x.  (46) 


Of  course  G*  is  unknown  so  assume  that  G  can  be  composed  into  an  ideal  portion  with  a  pertur¬ 
bation,  that  is  G  =  G*  +  AG.  The  system  input  then  becomes  u  =  G*y  +  AGy.  The  closed  loop 
system  then  appears  as, 

x  =  Acx  +  BAGy 

'T'  (47) 

y  =  Cx 

Here  a  Lyapunov  energy  function  can  be  defined  as  V(x)  =  \xT  Px  with  P  a  positive  definite, 
symmetric  matrix.  The  derivative  along  a  trajectory  is  given  by  V (x)  =  W/(x)  =  xTP  [Acx  +  5w]. 
With  some  algebraic  manipulation  the  following  relation  can  be  found, 


y(x)  =  l-xT  [PAc+A^P]x  +  xTPBxv. 


The  famous  Lyapunov  equation  is  seen  here 


PAc+AtcP=-Q, 


(48) 


(49) 


which  states  that  if  a  solution  (—0  to  the  matrix  equation  above  exists  then  the  derivative  along 
trajectories  will  be  less  than  zero  (  V (x)  <  0)  and  the  resulting  equilibrium  point  will  be  stable. 
Equation  48  can  be  expressed  as, 


V  (x)  =  —  ^xr<2x  +  yrw  with  PB  =  CT .  (50) 

This  energy  function  obviously  shows  a  dissipation  term  (—  Qx)  and  a  generation  term  (yrw). 
The  goal  of  our  adaptive  system  will  be  to  cancel  out  this  external  power  term  with  feedback 
control.  That  is  AG  needs  to  be  designed  such  that  the  V  (AG)  =  — yrw,  to  ensure  we  have  a  stable 
system. 

V (xAG)  =  V (x)  +  E (AG)  =  —  ^xr<2x  +  yrw  — yrw  <  0  V  x.  (51) 
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Define  a  new  energy  function  V(AG)  =  \tr(AGy  1  AGr)  with  y  >  0.  The  derivative  along 
trajectories  is  calculated  to  be 

V(AG)  =  tr(AG'y~1AGT).  (52) 

Now,  it  was  found  to  be  that  AG  =  — yy T  y  so  that, 

V  (AG)  =  tr(-yyTyy~1AGT)  =  -tr(-yyT  AGT)  =  -tr(  wTy)  =  -wry.  (53) 

Which  satisfies  our  condition  in  Equation  5 1 .  In  summary  adaptive  feedback  control  will  be  stable 
if  the  following  criteria  is  satisfied. 


PAc+AtcP  =  -Q 
PB  =  Ct 


(54) 


These  conditions  are  also  satisfied  if  the  system  is  strict  positive  real  (SPR)  with  no  non-minimum 
phase  zeroes.  Simply  written  as, 


A spr  <^=4>  CB  >  O&Minimum  phase  open  loop  system 


(55) 


This  process  is  done  exactly  the  same  for  discrete  time  systems;  the  only  difference  is  that  for 
strict  positive  realness,  the  so  called  Kalman- Yacubovic  equations  Fuentes  and  Balas  [2000]  must 
be  satisfied.  Which  are  (for  some  £  >  0), 

AtPA-P  =  -Q  =  -  (2 £P  +  LtL) 

AtPB  =  Ct  -LrW  (56) 

d+dt  =  WTW  +btpb. 


The  main  difference  here  is  that  the  D  term  cannot  be  zero  for  strict  positive  real  discrete  system. 
The  analysis  above  for  continuous  time  is  analogous  to  the  discrete  time  system. 
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5  Applications 

5.1  Cylinder  Wake 

The  DPOD  ANN-ARX  adaptive  control  design  approach  was  first  successfully  used  on  a  much 
simpler  flow  field,  the  cylinder  wake.  The  overall  goal  of  this  project  was  to  actuate  the  cylin¬ 
der  through  an  oscillating  displacement  to  reduce  the  Von-Karman  vortex  street.  Decreasing  the 
strength  of  these  vortices  reduces  the  amount  of  energy  put  into  the  wake  of  the  flow;  thus,  mini¬ 
mizing  the  drag  upon  the  cylinder.  This  project  was  taken  from  theoretical  formulation  to  numerical 
simulation  to  reduced  order  modeling  techniques  to  experimentation  validation.  The  cylinder  wake 
allowed  for  a  benchmark  flow  for  the  development  of  this  dynamic  model/control  design  approach. 
More  comprehensive  results  can  be  seen  in  Fagley  et  al.  [2009],  Siegel  et  al.  [2008], 

The  Von-Karman  vortex  street  creates  an  oscillating  lift  force  on  the  cylinder.  This  lift  force 
well  represents  the  natural  vortex  shedding  frequency  and  strength.  The  natural  shedding  frequency 
for  this  wake  was  approximately  5.6  Hz.  Figure  5  shows  a  forcing  simulation,  with  a  frequency 
equal  to  that  of  the  natural  shedding  frequency,  of  the  cylinder  wake  in  which  the  forcing  was 
started  at  exactly  1 80°  out  of  phase  from  the  lift  force.  As  seen  in  the  figure,  the  lift  force  goes 
through  a  transient  period  whilst  the  drag  decreases.  The  flow  then  begins  to  lock  in  with  the 
forcing  and  the  drag  returns  to  the  initial  value.  This  transient  phenomenon  is  important  to  see 
because  it  shows  that  open  loop  forcing  can  produce  a  desirable  flow  state  for  a  short  period  of 
time;  thus,  giving  feedback  control  a  promising  outcome. 


Figure  5:  Actuation  of  cylinder  wake  180°  out  of  phase.  The  resulting  transient  period  is  seen 
between  3.25  to  4s.  The  reduction  of  drag  is  directly  related  to  the  magnitude  of  the  limit  cycle  of 
the  lift  force. 

Forcing  parameter  cases  as  seen  in  Figure  8  were  simulated  with  CFD  software  such  as  Cobalt 
Solutions.  The  simulations  were  then  put  through  the  DPOD  process.  The  DPOD  spatial  modes 
are  seen  in  Figure  6  and  the  DPOD  time  coefficients  are  seen  in  Figure  7.  The  spatial  modes 
in  the  first  column  of  Figure  6  are  the  main  modes,  that  is  they  are  a  result  of  the  first  POD 
decomposition.  The  spatial  modes  in  the  second  column  of  Figure  6  are  the  shift  modes  which 
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Figure  6:  A  transient  forced  DPOD  spatial  mode  set  using  one  shift  mode  for  each  main  mode,  the 
first  3x2  DPOD  modes  are  shown.  Iso-contours  of  streamwise  velocity  are  shown,  solid  lines  are 
positive,  dashed  lines  negative. 


show  the  transient  affect  as  discussed  above.  The  mode  set  was  reduced  to  3  x  2  set  based  on 
energy  profile  analysis.  This  mode  set  contains  nearly  99.8%  of  the  flow  Siegel  et  al.  [2008].  Once 
the  DPOD  temporal  and  spatial  mode  sets  are  formulated,  they  are  validated  for  off  design  flow 
cases.  Errors  for  this  validation  are  on  the  order  of  0.5%  — >■  1%.  This  is  a  very  acceptable  range  of 
error  percentages.  Once  the  DPOD  mode  set  is  validated,  the  time  coefficients,  as  seen  in  Figure  7, 
for  the  selected  training  data  are  then  concatenated  and  the  dynamic  ANN-ARX  model  is  trained. 
Choosing  an  adequate  training  data  is  an  important  modeling  step  for  feedback  flow  control.  At 
first  training  data  in  which  the  flow  locked-in  with  forcing  actuation  was  used  to  formulate  the 
model.  This  model  had  poor  close  loop  dynamical  behavior  and  the  training  data  needed  to  be 
revised.  More  cases  were  added,  as  seen  in  Figure  8  and,  more  importantly,  more  transient  cases 
were  added  in  which  the  flow  field  didn’t  lock  in  with  the  actuation.  These  lower  amplitude  forcing 
cases  produced  very  long,  non-linear  transients.  The  model  then  accurately  predicted  closed  loop 
behavior  as  discussed  later. 

The  neural  network  topology  is  an  important  factor  for  realizing  the  non-linear  behavior  of  the 
data  to  be  identified.  Typically,  neural  networks  consist  of  three  layers:  input,  hidden  and  output. 
Activation  functions  are  also  determined  for  each  of  the  neurons.  The  sigmoid  function  is  the  most 
commonly  used  function  due  to  the  continuity  of  the  function  for  back  propagation  derivation  train¬ 
ing  algorithms.  Other  parameters  are  also  important  for  the  functionality  of  the  network  model. 
The  number  of  neurons  allows  the  network  to  fit  more  complicated  non-linear  trends.  Eventhough, 
the  more  neurons  a  network  contains  the  more  likely  the  model  will  be  overtrained  and  poorly 
simulate  off-design  data.  Also,  training  times  will  be  greatly  increased.  The  number  of  past  inputs 
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i  =  1,  j  =  1  i  =  1,  j  =  2 


15  20  25  30  35  40  45  15  20  25  30  35  40  45 

t/T  t/T 


Figure  7:  Mode  Amplitudes  of  the  open  loop  forced  simulation  j^  =  1  and  ^  =  0.25.  Forcing 
activated  at  ^  =  18  and  stopped  at  ^  =  33,  after  15  full  forcing  cycles.The  first  3x2  DPOD  mode 
amplitudes  are  shown. 


Figure  8:  Training  data  set  used  for  ANN-ARX  identification.  Lower  amplitude  forcing  near 
natural  frequency  produced  slow,  large  transients. 
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Table  1:  ANN-ARX  3x2  model  parameters 


Input 

#  Past  Inputs 

Sampling  Delay 

Total  Time  History 

Re  Number 

1 

10 

10 

Actuator  Position 

4 

2 

8 

flt.i 

1 

1 

1 

«2,1 

3 

8 

24 

<33,1 

3 

8 

24 

<31,2 

1 

1 

1 

<32,2 

1 

12 

12 

03,2 

1 

12 

12 

and  time  histories  of  previous  simulated  outputs  also  is  a  crucial  design  factor.  A  new  feature  was 
added  which  allows  for  a  time  tapped  delay  so  every  nth  point  will  be  sampled.  This  gives  for 
a  much  larger  time  history  while  keeping  the  number  of  inputs  lower  which  yields  much  shorter 
training  times.  All  of  these  parameters  were  adjusted,  trained,  and  repeated  until  a  proper  model 
was  found.  The  parameters  can  be  seen  in  the  Table  1. 


=  1.1-1 


1-1. 1-2 


Figure  9:  Off  design  validation  of  ANN-ARX  model.  Forcing  case  for  =  0  and  jj~!  =  10%. 
Red  lines  show  the  simulation  result  of  the  model  and  the  blue  lines  are  the  actual  CFD  data. 

The  model  is  also  validated  for  off  design  actuation.  As  shown  below  in  Figure  9,  the  3x2  ANN- 
ARX  model  does  an  excellent  job  of  simulating  the  mode  amplitudes  for  an  off  design  actuation 
case  with  forcing  of  -jr2-  =  0  and  =  10%. 

The  adaptive  feedback  control  algorithm  in  section  2.2.1  was  designed  and  mode  ci2,\  was 
fed  back  to  the  vertical  cylinder  position.  Once  the  gain  aggressiveness  (y)  in  equation  2.25  was 
properly  conditioned  and  correspondingly  good  results  were  seen  with  the  ANN-ARX  closed  loop 
simulation  (i.e.  reduction  of  mode  ci2.  \ )  were  seen  the  same  control  algorithm  was  plugged  into 
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a  CFD  simulation.  The  CFD  simulation  shows  qualitatively  accurate  results  when  compared  to 
the  ANN-ARX  model.  Thus  proving  that  the  model  does  capture  the  close  loop  dynamics  of  the 
cylinder  wake  flow  field. 


Figure  10:  Closed  loop  simulation  of  the  3x2  ANN-ARX  model(left).  CFD  closed  loop  simula- 
tion(right).  Similar  dynamics  are  shown  between  each  of  the  models. 

The  CFD  simulation  allows  for  the  ability  to  analyze  flow  characteristics  such  as  pressures, 
velocities,  densities,  etc.  Here  the  surface  pressures  of  the  cylinder  can  be  integrated  around  the 
surface  to  give  resulting  lift  forces  and  more  importantly  drag  force.  Figure  1 1  shows  a  reduction 
of  drag  on  the  cylinder  up  to  16%.  Thus,  proving  the  control  development  strategy  successful  for 
simple  flow  fields. 
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Drag  on  the  Cylinder  Fx  Cylinder  Position  /  Controller  Output 


Figure  11:  CFD  closed  loop  simulation  using  direct  adaptive  control.  Drag  force  reduction  (left) 
and  actuation  input/cylinder  position  (right). 
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5.2  Aero  Optics 

Aero-optical  systems  are  designed  for  the  transmission  of  light  beams  through  aerodynamic  flow 
fields.  For  optical  design  purposes,  the  flow  has  to  be  considered  as  a  time-varying  optical  ele¬ 
ment  due  to  density  variations  of  the  fluid  in  the  beam  path,  which  lead  to  changes  in  the  optical 
path  length  (OPL).  The  correlation  between  density  fluctuations  and  variations  in  the  gas  refrac¬ 
tive  index,  whose  integral  is  the  optical  path  length,  is  described  by  the  Gladstone-Dale  relation, 
which  states  a  linear  relationship  where  the  proportionality  constant  is  the  Gladstone-Dale  constant 
McMackin  et  al.  [1997],  Jumper  and  Fitzgerald  [2001], 

Aero-optical  aberrations  can  be  loosely  grouped  into  two  categories.  The  first,  associated  with 
the  large  scales  in  the  flow  field,  includes  boresight  (tracking)  errors.  Some  of  these  errors  could 
be  alleviated  by  current  adaptive  optical  systems  because  of  their  relatively  large  length  scales  and 
slow  time  scales.  The  second  category  includes  errors  such  as  beam  spreading,  scintillation  and 
reduction  of  resolution,  contrast,  etc.,  due  to  the  small  scale  turbulent  motion  in  the  flow.  Although 
recently  there  have  been  some  indications  that  the  large  scale  motion  can  also  cause  errors  typically 
associated  with  small  scales,  this  classification  still  seems  appropriate.  In  the  past,  corrections  for 
these  aberrations  in  aero-optical  systems  were  based  on  compensation  using  optical  components. 
This  approach  led  to  the  development  of  highly  complex  adaptive  optical  systems.  While  truly 
impressive  results  have  been  achieved  for  telescopes  used  in  astronomy,  which  have  to  correct  for 
aberrations  due  to  the  earth’s  atmosphere,  state-of-the-art  adaptive  optical  systems  have  only  met 
limited  success  on  airborne  applications.  To  a  large  extent,  this  is  due  to  the  vastly  different  length 
and  time  scales,  compared  to  terrestrial  astronomy  applications,  present  on  envisioned  airborne 
optical  platforms.  One  of  the  problems  is  that  all  adaptive  optical  systems  rely  on  mechanically 
moving  some  component  (usually  a  mirror  surface)  to  adapt  to  beam  distortions  due  to  the  flow 
field.  Current  systems  are  still  limited  in  their  bandwidth,  frequency,  and  field-of-vision  and  are 
unable  to  correct  for  disturbances  of  all  length  and  time  scales.  The  same  limitations  apply  to  wave 
front  sensors,  and  although  new  wave  front  sensors  are  being  developed  to  analyze  distortions  in 
aero-optical  applications,  they  still  suffer  from  limitations  that  make  them  only  marginally  usable 
in  environments  comparable  to  the  one  in  airborne  applications  Trolinger  et  al.  [2005]. 

A  relatively  new  approach  to  understanding  and  controlling  the  optical  aberrations  observed  in 
applications  is  to  look  at  directly  controlling  the  flow  field  to  minimize  the  strength  of  the  structures 
responsible  for  the  optical  distortions.  If  control  is  possible,  and  recent  research  Seidel  et  al. 
[2005],  Siegel  et  al.  [2005a]  has  shown  initial  successes  in  this  research  area,  the  aberrations  could 
be  reduced,  possibly  overcoming  a  significant  hurdle  on  the  way  to  the  implementation  of  airborne 
optical  systems.  In  this  context,  two  fundamentally  different  ideas  of  flow  control  are  currently 
being  explored.  The  first  is  open-loop  active  flow  control  (AFC),  which  introduces  small  amplitude 
disturbances  at  sensitive  locations  in  the  flow,  which  will  then  be  amplified  using  the  flow’s  inherent 
instabilities  and  thus  can  lead  to  large  global  changes  in  flow  behavior.  The  best  known  example 
of  this  technique  is  the  delay  of  separation  on  airfoils  at  large  angles  of  attack  by  introducing 
small  disturbances  upstream  of  the  separation  point  that  keep  the  flow  attached  by  energizing 
the  boundary  layer.  However,  even  AFC  is  not  able  to  react  dynamically  to  changing  operating 
conditions.  Furthermore,  AFC  relies  on  creating  disturbances,  which  in  the  context  of  aero-optics 
applications  is  detrimental  to  improving  system  performance.  The  problems  outlined  above  led  to 
the  development  of  a  second  idea  for  flow  control:  feedback  (closed-loop)  flow  control.  The  key 
components  of  the  feedback  system  are  flow  sensors,  the  state  estimator,  a  controller  module  and 
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the  actuators.  The  information  obtained  from  sensors,  which  observe  instantaneous,  localized  flow 
quantities,  is  analyzed  in  the  state  estimator.  The  reduced  data  is  then  supplied  to  the  controller 
module  to  determine  an  output  signal  that  is  used  to  drive  the  actuators.  It  has  been  shown  that 
using  feedback  flow  control,  it  is  possible  to  achieve  goals  such  as  the  suppression  of  the  von 
Karman  vortex  street  Siegel  et  al.  [2005a]  that  are  not  achievable  using  active  open  loop  control. 
Also,  feedback  flow  control  is  typically  more  efficient  than  active  flow  control  due  to  its  ability  to 
optimize  the  forcing  input  to  match  the  current  flow  state.  Furthermore,  feedback  flow  control  is 
tolerant  against  environmental  changes,  since  it  senses  the  actual,  instantaneous  flow  field  rather 
than  operating  on  assumed  states.  And  finally,  because  feedback  control  is  used  on  demand,  it  does 
not  have  a  detrimental  effect  in  other  flow  regimes.  These  properties  make  feedback  flow  control 
superior  to  active  flow  control  despite  the  increased  complexity. 

For  the  aero-optical  problem,  two  different  approaches  for  implementing  AFC  have  been  in¬ 
vestigated  to  date.  The  first,  developed  by  Jumper  and  co-workers  [see  e.g.  Gordeyev  et  al.,  2005], 
is  based  on  the  idea  that  regularizing  the  flow  will  yield  a  flow  field  that  is  more  deterministic 
with  respect  to  its  unsteadiness.  This  is  achieved  by  forcing  the  flow  with  a  known  disturbance 
signal  (frequency,  phase,  and  amplitude).  The  knowledge  of  the  resulting  structure  of  the  flow 
field,  including  the  (approximate)  strength  and  phase  of  the  large,  most  optically  active  struc¬ 
tures,  simplifies  the  task  of  the  adaptive  optics  system,  which  targets  the  now  known  distortions. 
While  this  approach  has  been  shown  to  yield  good  results,  it  also  strengthens  the  large  structures, 
which  is  counterproductive  when  the  goal  is  to  minimize  their  optical  aberrations.  The  second 
approach,  which  was  developed  by  Glezer  and  co-workers  Oljaca  and  Glezer  [1997],  Vukasinovic 
et  al.  [2004],  hinges  on  the  observation  that  the  large,  and  therefore  low  frequency,  coherent  struc¬ 
tures  in  the  flow  field  can  be  destroyed  by  high  frequency  forcing.  Although  they  have  shown 
successful  reduction  of  the  large  structures  and  their  associated  optical  distortions,  it  stands  to 
reason  that  increasing  the  energy  contained  in  the  small  structures  is  detrimental  to  the  optical 
performance  because  it  strengthens  exactly  those  structures  that  provide  the  high  frequency,  small 
scale  aberrations  that  are  outside  the  realm  of  correction  of  current  adaptive  optical  systems. 

In  contrast  to  the  above  mentioned  open-loop  AFC  research,  feedback  flow  control  was  used  in 
this  project  to  control  the  unsteady  structures  in  the  flow  field.  The  underlying  assumption  is  that  if 
the  large,  coherent  structures  can  be  successfully  weakened,  their  aberrations  will  diminish  as  well. 
In  addition,  when  energy  is  extracted  from  the  large  structures,  the  energy  available  to  create  small 
scales  is  diminished,  delaying  the  development  of  turbulence.  This  notion  has  been  successfully 
demonstrated  for  the  Karman  vortex  street,  where  successful  control  of  the  shedding  frequency 
also  yielded  an  amplitude  reduction  in  higher  harmonics  Siegel  et  al.  [2003a,  2005a],  It  is  this 
combination  of  effects  that  holds  the  promise  of  successfully  controlling  the  optical  aberration  due 
to  the  flow  over  the  aperture  of  airborne  optical  platforms. 

As  outlined  above,  systems  currently  in  use  suffer  from  aberrations  that  are  outside  the  capabil¬ 
ities  of  state-of-the-art  adaptive  optics  systems.  To  alleviate  the  tasks  to  be  borne  by  the  adaptive 
optical  system  has  the  potential  to  provide  a  highly  sought  after  functionality.  The  renewed  interest 
in  these  systems,  with  much  more  stringent  requirements  for  accuracy,  requires  novel  methods  to 
reduce  the  detrimental  effect  of  the  flow  over  the  aperture. 
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Figure  12:  Turbulent  shear  layer  Van  Dyke  [1982]. 


5.2.1  Technical  approach 

The  core  technologies  for  this  aero-servo-optics  project,  aimed  at  reducing  optical  aberrations  in 
airborne  platforms,  were  derived  from  a  combination  of  the  fluid  dynamics,  controls,  and  optics 
research  areas,  with  a  strong  emphasis  on  fluid  dynamics  and  controls.  The  optical  performance  of 
the  proposed  aerodynamic  control  system  was  judged  using  the  well-known  correlation  between 
the  fluid  density  and  index-of-refraction  fields  in  gases.  The  underlying  assumption  was  that  as 
long  as  reliable  data  of  the  density  field  is  available  (through  simulations  or  experiments),  a  good 
representation  of  the  optical  properties  should  be  achievable. 

From  a  fluid  dynamics  point  of  view,  shear  layers  develop  large,  coherent  structures  due  to 
their  inherent  inviscid  instability.  As  shown  in  Figure  12,  there  are  several  stages  in  this  devel¬ 
opment.  First,  large,  laminar,  S-shaped  structures  are  generated  behind  a  splitter  plate.  After  the 
initial  growth  of  these  structures,  smaller  structures  develop  due  to  further  instability  mechanisms, 
eventually  leading  to  small  scale  turbulence.  However,  as  shown  in  the  figure,  even  when  the  flow 
is  turbulent,  the  large  coherent  structures  persist.  It  is  these  structures  that  are  responsible  for 
the  large  boresight  errors  in  aero-optical  applications.  However,  because  of  their  large  size  and 
relatively  low  frequency,  an  adaptive  optical  system  can  correct  for  their  effect  if  their  amplitude 
and  phase  are  known.  Another  important  aspect  to  note  in  Figure  12  is  that  the  turbulent  motion 
develops  as  a  consequence  of  the  primary  flow  instability  that  creates  the  S-shaped  shear  layer 
structures.  This  is  in  accordance  with  observations  of  the  energy  cascading  from  large  to  small 
scales  Pope  [2000], 

For  the  aero-servo-optics  project,  this  energy  cascade  is  of  crucial  importance:  If  feedback 
flow  control  is  able  to  reduce  the  strength  of  the  primary  shear  layer  structures,  it  follows  that 
there  is  less  energy  in  these  disturbances  and  therefore  less  energy  is  available  to  generate  small 
scale  turbulent  structures,  resulting  in  significantly  reduced  optical  aberrations  on  all  scales.  An 
effect  similar  to  the  one  outlined  above,  namely  the  reduction  of  the  amplitude  of  high  frequency 
disturbances  by  controlling  a  lower  frequency,  has  been  observed  in  low  Reynolds  number  wake 
flows  Cohen  et  al.  [2003a]. 

There  were  certain  aspects  of  the  investigations  that  lent  themselves  to  investigations  using 
simulations,  while  other  aspects  were  explored  more  easily  and  efficiently  in  the  experiment.  For 
example,  only  small  portions  of  the  flow  can  be  measured  simultaneously  in  the  experiment  since 
the  sensing  options  are  limited.  Thus,  meaningful  sensor  locations  can  be  derived  much  more 
easily  from  simulation  results  since  the  entire  flow  field  data  with  all  its  variables  is  available. 
However,  in  order  to  vary  flow  and  actuation  parameters,  a  new  simulation  needs  to  be  performed 
for  each  set  of  parameters,  with  the  associated  cost  and  time  requirements.  In  contrast,  once  the 
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experimental  hardware  is  in  place,  flow  and  actuation  parameters  can  be  varied  easily  and  quickly. 
Thus  this  type  of  investigation  is  best  performed  using  experiments.  In  summary,  by  using  both 
experiment  and  computation  in  parallel,  the  fastest  possible  progress  was  achieved. 

To  provide  an  overview  of  the  feedback  flow  control  design  cycle  used  in  this  research  project. 
Figure  1  shows  the  main  building  blocks  in  the  process.  The  development  started  with  building  a 
database  of  flow  states  based  on  CFD  simulation  results.  First,  the  natural  (i.e.  without  any  control 
input)  flow  field  was  simulated.  Then,  a  number  of  simulations  were  performed  where  periodic 
blowing  and  suction  was  used  to  introduce  disturbances  at  a  given  frequency  and  amplitude  into 
the  flow  (see  Section  5.2. 2.4).  The  results  of  all  these  simulations  were  analyzed  using  Proper 
Orthogonal  Decomposition  (POD),  which  resulted  in  POD  spatial  modes  as  well  as  the  POD  time 
coefficients  for  each  time  step  of  all  simulations. 

These  POD  modes  and  time  coefficients  were  then  used  for  the  development  of  a  reduced  order 
model  (ROM).  In  the  present  effort,  a  wavenet  ARX  topology  was  chosen  (see  Section  4.2.2. 1). 
Once  the  model  performance  was  validated  against  the  original  CFD  data,  state-of-the-art  feedback 
controller  design  tools  were  used  to  develop  a  controller.  Iterative  testing  (iteration  loop  “1”  in 
Figure  1)  lead  to  a  controller  design  that  achieved  the  predefined  control  goal  of  minimizing  the 
optical  distortion  for  a  given  aperture. 

In  addition,  the  POD  spatial  modes  were  scrutinized  for  flow  state  estimation  purposes.  Sensors 
placement  studies,  which  were  used  to  determine  the  number  and  locations  of  flow  sensors,  were 
performed  using  the  computed  flow  quantities  on  the  wall  behind  the  backward  facing  step  (see 
Section  4.2.4.2).  With  the  sensors  chosen,  a  flow  state  estimator  was  developed.  This  estimator 
determined  the  global  flow  state  based  on  the  sensor  readings,  i.e.  it  established  field  data  from 
only  the  sensor  information. 

At  this  stage,  the  flow  state  estimator  and  the  controller  were  introduced  into  the  CFD  simu¬ 
lations  and  feedback  controlled  simulations  were  performed.  The  results  of  this  simulation  were 
scrutinized  to  investigate  the  effect  of  the  control  input  on  the  flow  field,  as  well  as  their  effect  on 
the  overall  figure  of  merit,  the  OPD.  As  indicated  in  Figure  1,  multiple  iteration  paths  were  open  at 
this  point.  Path  “2”  in  Figure  1  could  be  taken  if  the  results  indicated  that  the  controller,  designed 
using  the  wavenet  model,  is  not  performing  well.  It  could  also  be  possible  that  the  findings  indi¬ 
cate  problems  with  the  wavenet  itself,  which  would  be  remedied  using  iteration  path  “3”.  Finally, 
the  research  design  is  flexible  enough  to  also  allow  inclusion  of  feedback  controlled  data  into  the 
POD  database  to  improve  the  fidelity  of  both  model  and  controller  development.  All  steps  outlined 
above  will  be  described  in  detail  in  this  report. 

5.2.1. 1  Basic  flow  parameters  Several  parameters  had  to  be  considered  to  arrive  at  flow  con¬ 
ditions  that  on  the  one  hand  result  in  large  enough  (i.e.  measurable)  density  changes,  but  that,  on 
the  other  hand,  do  not  result  in  vortex  shedding  frequencies  that  were  too  high  to  be  measured  suc¬ 
cessfully.  In  addition,  the  slower  the  flow,  the  easier  it  is  to  implement  a  feedback  control  system. 
These  considerations  resulted  in  the  following  key  experimental  parameters: 

•  Mach  number  Ma  =  0.3  at  the  inflow 

•  Step  height  H  =  0.15m 

•  Free  stream  velocity  at  step  [/«,  =  140m/.v 
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5.2.1.2  Optical  definitions  The  optical  properties  of  the  fluid  forming  the  shear  layer  behind 
the  backward  facing  step  were  evaluated  based  on  the  Gladstone  Dale  relation 

n(x,y,z,t)  =  l+kGDp(x,y,z,t),  (57) 

where  n  is  the  fluid’s  index  of  refraction,  kc,D  is  the  Gladstone-Dale  constant,  kc,D  =  2.289  x 
10_4m3/kg,  and  p  is  the  fluid  density.  The  optical  path  length  (OPL)  can  be  obtained  by  integrating 
the  index  of  refraction  along  the  beam  path,  L, 

OPL(x,y,z,t)  =  f  n(x,y,z,t)  dl  =  L  +  kcD  [  p(x,y,z,t)  dl.  (58) 

Jo  Jo 

Since  the  differences  in  OPL  over  a  given  aperture  are  typically  on  the  order  of  the  wave  length  of 
the  beam,  it  is  common  to  express  the  wave  front  distortion  as  the  optical  path  difference  (OPD), 
which  is  defined  as  the  local,  instantaneous  OPD  minus  the  spatial  mean  over  the  aperture.  As¬ 
suming  the  beam  propagates  in  the  y-direction,  this  can  be  expressed  as 

OPD(v,z,t;y)  =  OPL(x,z,t;y) -OPL(x,z,t;y)xz.  (59) 

These  equations  were  used  when  analyzing  the  CFD  results  because  the  density  field  p(x,y,z,t)  is 
computed  directly  from  the  governing  equations. 

5.2.2  Numerical  Simulation 

The  exact  geometry  for  the  simulations  of  a  free  shear  layer  was  developed  in  conjunction  with  the 
design  of  the  experiment.  Comparisons  to  experimental  data  performed  at  USAFA  aeronautics  lab¬ 
oratory  of  the  unforced  flow  field  data  were  performed  to  validate  the  accuracy  of  the  simulations 
and  to  judge  the  necessary  grid  resolution  to  resolve  the  relevant  flow  features.  In  particular,  the 
optical  path  difference  (OPD)  was  used  as  the  main  optical  figure-of-merit.  It  should  be  noted  that 
once  a  time  accurate,  spatial  density  distribution  is  available  from  the  computations,  calculating 
the  index  of  refraction  field  and  the  resulting  OPD  is  possible  with  a  small  computational  effort 
compared  to  the  CFD  simulations.  During  the  course  of  this  study,  Mani  et  al.  [2008]  published 
an  article  outlining  the  resolution  requirements  determined  by  their  aero-optical  simulations.  They 
concluded  that  the  resolution  requirements  for  an  aero-optic  simulation  match  the  ones  for  a  well 
resolved  LES  simulation. 

To  build  a  database  of  flow  states  that  would  be  used  to  define  the  reduced  order  model  for 
the  flow  field,  unforced  simulations  were  performed  first.  In  a  second  step,  open  loop  active  flow 
control  (AFC),  which  in  the  simulations  was  implemented  using  an  externally  controlled  blowing/- 
suction  boundary  condition  (see  below),  was  studied  and  the  data  was  added  to  the  development 
cycle  of  the  database.  These  forcing  cases  were  particularly  valuable  for  describing  the  transient 
flow  features  present  during  the  initial  development  of  the  open  loop  forced  shear  layer  as  well 
as  the  vortex  pairing  that  occurred  when  forcing  was  initiated.  The  results  from  the  simulations 
provided  a  comprehensive  database  of  the  free  shear  layer,  which  was  used  to  develop  feedback 
control  strategies  as  well  as  to  compare  the  effectiveness  of  feedback  control  applied  to  the  aero- 
optics  problem. 
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(a)  2D  Grid 


(b)  B/S  slot 


Figure  13:  a)Two-dimensional  CFD  grid,  b)  Grid  at  the  step  showing  the  blowing/suction  slot. 


5.2.2.1  Grid  generation  The  experimental  ramp  geometry  was  used  in  the  simulations.  The 
grids  for  the  simulations  were  generated  using  the  SimCenter  software  developed  at  Mississippi 
State  University  Marcum  and  Weatherill  [1995],  Since  the  geometry  is  essentially  two-dimensional, 
a  planar  grid  was  generated  first  and  then  this  grid  was  extruded  in  the  spanwise  direction.  The 
step  height  was  H  =  0. 15m,  and  the  ramp  length  Lr  =  0.85m.  To  ensure  that  no  disturbances  reach 
the  outflow  boundary,  the  domain  length  was  Lx  =  4m  downstream  of  the  step.  The  domain  height 
matched  the  experimental  setup  in  the  USAFA  wind  tunnel,  Ly  =  0.85m. 

The  main  difference  to  the  experimental  geometry  was  that  the  forcing  chamber  (see  Figure  19) 
was  not  included  in  the  simulations.  Instead,  only  a  short  section  of  the  slot  ( b  =  1mm)  was  mod¬ 
eled  and  a  blowing  and  suction  boundary  condition  was  applied  at  the  base  of  the  slot.  Figure  13 
shows  the  final  two-dimensional  grid.  A  zoomed  view  of  the  step  with  the  blowing/suction  slot  is 
given  in  Figure  13.  The  grid  spacing  at  the  step  was  defined  to  be  Ax  =  0.1mm.  This  grid  contains 
approximately  58,000  nodes  and  90,000  elements.  Grid  clustering  was  used  on  the  bottom  wall 
and  in  the  region  of  interest  in  the  free  shear  layer.  In  order  to  resolve  the  blowing  and  suction  slot 
geometry,  the  grid  also  had  to  be  refined  near  the  step  edge.  The  boundary  layer  grid  spacing  was 
chosen  such  that  the  final  y+  value  at  the  step  was  y+  ~  1. 

For  the  three-dimensional  simulations,  this  grid  was  extruded  in  the  spanwise  direction.  The 
spanwise  step  size  was  chosen  as  A z  =  1mm.  Solutions  at  various  domain  widths  (L-/H  =  1 , 2, 3, 4) 
were  computed  to  ensure  that  the  pertinent  shear  layer  dynamics  were  captured  in  the  simulations. 
Figure  14  shows  the  geometry  for  the  case  Lz/H  =  2. 

In  addition  to  the  grid  for  the  CFD  simulations,  a  grid  for  the  beam  propagation  was  devel¬ 
oped.  The  beam  domain  size  was  chosen  as  0  <  x/H  <  3,  —1  <  y/H  <  0.5,  and  spanned  the 
whole  spanwise  domain,  which  represents  a  sufficiently  large  domain  to  investigate  various  optical 
apertures  while  being  able  to  maintain  reasonable  resolution.  46  x  81  x  41  points  were  used  in 
the  x-,  y-,  and  z-directions,  respectively.  The  grid  was  designed  as  a  structured  grid  with  one  grid 
direction  aligned  with  the  predominant  beam  direction.  This  approach  facilitated  the  computation 
of  the  OPL  and  OPD  (see  Equations  58  and  59)  since  the  path  integral  was  along  grid  lines.  For 
the  interpolation  of  the  CFD  data  onto  the  beam  grid,  the  ’tap’  capabilities  in  COBALT  were  used. 
Taps  were  initially  designed  as  measurement  locations,  but  for  the  current  research,  using  taps  to 
extract  the  flow  field  on  the  beam  grid  ensured  that  the  numerical  methods  for  integration  of  the 
Navier-Stokes  equations  and  interpolation  were  consistent. 
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Figure  14:  Three  dimensional  domain  with  CFD  grid,  Lz/H  =  2. 

5. 2.2.2  Results  Figure  ??  shows  the  comparison  of  the  mean  flow  u-velocity  profiles  at  x/H  = 
0, 0.5, 1 , 2, 3, 4  obtained  from  experiments  (symbols)  and  simulations  (lines).  The  separating  bound¬ 
ary  layer  had  a  thickness  of  899  ~  8mm,  which  for  a  fully  turbulent  boundary  layer  corresponds 
to  Reg  ~  8500.  As  the  shear  layer  develops,  at  downstream  positions  x/H  <  2,  the  shear  layer  in 
the  simulations  is  not  spreading  as  quickly  as  measured  in  the  experiments.  Further  downstream 
at  x/H  =  3  the  profiles  are  in  very  good  agreement  and  at  x/H  =  4  the  simulation  results  show 
a  slightly  larger  shear  layer  thickness.  This  increased  spreading  rate  in  the  simulation  data  was 
attributed  to  the  behavior  of  the  DES  turbulence  model  as  the  separating  flow  transitions  from  a 
RANS  based  boundary  layer  calculation  to  a  Large  Eddy  Simulation.  The  initial  lack  of  structures 
in  the  flow  led  to  reduced  shear  layer  growth.  As  the  structures  developed,  the  growth  rate  matched 
the  experimental  and  theoretical  results  well,  but  further  downstream,  the  grid  resolution  is  insuffi¬ 
cient  to  maintain  the  coherent  structures  in  the  flow  and  numerical  diffusion  results  in  an  increased 
spreading  rate. 

Instantaneous  results  of  the  simulations  are  show  in  Figure  15.  In  Figure  15a,  the  flow  struc¬ 
tures  are  visualized  using  an  iso-surface  of  the  Q  vortex  identification  criterion  Jeong  and  Hussain 
[1995]  colored  by  pressure.  At  this  instant,  the  shear  layer  (Kelvin-Helmholtz)  vortices  are  starting 
to  form  approximately  one  step  height  downstream  of  the  separation  point,  with  increased  span- 
wise  coherence  one  wavelength  further  downstream.  The  instantaneous  isosurface  of  density  is 
shown  in  Figure  15b.  Comparing  the  Q-vortex  structures  with  the  density  isosurface  shows  that 
there  is  a  very  strong  correlation. 

Furthermore,  it  is  interesting  to  note  that  the  density  isosurface  shows  only  the  large  scale  struc¬ 
tures  while  suppressing  the  smaller  scales  at  the  step  as  well  as  in  the  recirculation  region  below 
the  shear  layer.  This  is  due  to  the  deeper  “pressure  well”,  and  the  concomitant  drop  in  density, 
inside  the  largest  structures.  For  feedback  flow  control  that  targets  the  coherent  motion  in  the  shear 
layer,  this  behavior  of  the  density  field  is  highly  desirable  because  density  behaves  like  a  filter 
and  density  isosurfaces  identify  the  flow  structures  of  interest.  In  addition,  since  the  optical  path 
length  is  a  linear  function  of  density,  density  is  in  fact  the  quantity  of  interest  for  the  aero-optics 
problem.  This  is  shown  in  Figure  16,  where  the  flow  structures,  identified  using  the  Q-criterion 
Jeong  and  Hussain  [1995],  are  shown  in  grey  and  the  OPD  is  plotted  in  color  at  the  top  of  the  beam 
grid.  Comparing  the  OPD  results  and  the  density  isosurface  plotted  in  Figure  15  shows  a  strong 
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(b)  p  =  0.95kg/m3  isosurface. 


(a)  Q  =  2  •  106  1/s2  isosurface. 


Figure  15:  Instantaneous  representation  of  the  structures  in  the  flow  field.  Iso  surfaces  colored  by 
pressure. 


Figure  16:  Instantaneous  flow  structures  visualized  using  an  isosurface  Q  =  5  •  105  1/s2  and  OPD 
(color). 

correlation  between  the  “valleys”  in  the  wave  front  (negative  values  of  OPD  shown  in  blue)  and  the 
location  of  the  flow  structures.  The  initial  vortex  shedding  frequency  for  the  shear  layer  can  be  esti¬ 
mated  from  theory  based  on  St  =  Fd/Uoo  =  0.012  Hasan  [1992]  to  be  Fn  ~  2000Hz.  Using  probes 
at  y/H  =  0  and  various  streamwise  positions,  the  frequencies  with  the  highest  amplitudes  ranged 
around  F  ~  400Hz  (Figure  17),  indicating  that  vortex  pairings  had  occurred  upstream  of  the  probe 
locations  [see  Seidel  et  al.,  2009,  for  more  details].  While  the  total  simulation  time  was  too  short 
for  a  detailed  spectral  analysis,  the  results  provided  a  good  indication  of  the  frequency  of  the  natu¬ 
rally  occurring  structures.  The  results  were  also  in  good  agreement  with  the  results  obtained  from 
the  experiments.  Because  the  extent  in  the  streamwise  direction  spanned  by  the  density  probes  was 
commensurate  with  the  area  of  interest  for  optical  performance,  this  frequency  provided  an  approx¬ 
imate  target  frequency  to  investigate  the  effect  of  open  loop  forcing  on  the  optically  relevant  shear 
layer  structures.  To  further  analyze  the  flow  field  and  the  structures  in  the  shear  layer,  the  simu¬ 
lation  data  was  reduced  using  Proper  Orthogonal  Decomposition  Sirovich  [1987],  Berkooz  et  al. 
[1993],  Holmes  et  al.  [1996].  The  u- velocity,  v- velocity,  and  density  were  analyzed  to  examine 
which  of  these  quantities  provided  a  meaningful  representation  of  the  flow  structures  with  a  focus 
on  the  coherent  structures  in  the  shear  layer.  Figure  18  shows  the  first  six  spatial  modes  (plotted 
in  pairs  to  show  the  traveling  vortex  nature  of  the  shear  layer  structures)  and  their  corresponding 
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Figure  17:  Density  spectra  between  x/H  =  0  and  x/H  =  2.6.  Unforced. 


time  coefficients  for  the  full  length  of  time  available  in  the  data  set.  Figure  18a  shows  that  the 
u-velocity  POD  modes  do  not  have  any  discernable  spatial  structure,  which  is  also  reflected  in  the 
lack  of  a  dominant  frequency  in  the  time  coefficients  (Figure  18b).  In  contrast,  the  first  two  pairs 
of  POD  modes  of  the  v-velocity  (Figure  18c  and  d)  exhibit  a  distinct,  highly  coherent  structure  that 
is  indicative  of  traveling  waves.  The  periodic  character  of  the  time  coefficients  corroborates  that 
these  POD  modes  capture  the  vortex  street  in  the  shear  layer.  The  same  holds  for  the  density  POD 
modes  (Figure  18e  and  f)  that  identify  the  shear  layer  structures  and  their  spanwise  distortions.  It 
is  interesting  to  note  that  the  POD  of  both  the  v-velocity  and  density  show  the  strongest  peaks  in 
the  shear  layer  while  the  smaller  scales  in  the  recirculation  zone  below  are  suppressed;  in  contrast, 
the  u-velocity  POD  modes  include  these  structures.  Furthermore,  the  dominant  mode  pairs  of  both 
the  v-velocity  and  density,  modes  1-2  and  3-4,  show  the  same  frequency  content  but  are  not  pe¬ 
riodic,  corroborating  the  fact  that  the  spectral  peak  around  the  dominant  frequency  is  very  broad 
in  the  natural  flow.  Another  important  point  to  note  is  the  strong  spanwise  coherence  of  modes 
1-2  of  both  the  v-velocity  and  density.  This  indicates  that  the  dominant  structures,  as  measured 
by  the  magnitude  of  the  singular  values  obtained  from  the  POD  procedure  (not  shown),  are  indeed 
spanwise  “rollers”,  even  in  the  unforced,  three-dimensional  flow.  Forcing,  as  described  in  the  next 
section,  only  increases  this  spanwise  coherence. 
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Table  2:  Summary  of  computed  forcing  cases. 


400Hz 

600Hz 

800Hz 

1000Hz 

A/Uoo  =  0.3 

X 

X 

X 

X 

A/Uoo  =  0.2 

X 

X 

X 

X 

A/Uoo  =  0.1 

X 

X 

X 

X 

A/Uoo  =  0.05 

X 

X 

X 

X 

A/Uoo  =  0.01 

X 

X 

X 

X 

While  this  POD  analysis  was  not  exhaustive  (other  quantities,  e.g.  vorticity,  could  be  analyzed), 
it  showed  that  either  the  v-velocity  or  the  density  are  well  suited  for  developing  a  reduced  order 
model  for  control  purposes  for  this  research  program.  From  an  aero-optics  perspective,  the  den¬ 
sity  is  clearly  the  best  quantity  for  model  and  feedback  controller  development  due  to  its  direct 
influence  on  the  optical  properties  of  the  flow. 

5.2.2.3  Actuation  Since  the  forcing  chamber  was  not  part  of  the  backward  facing  step  flow 
geometry,  it  was  modeled  in  a  separate  simulation  to  verify  that  the  slot  exit  velocity  was  uniform 
in  the  spanwise  direction.  The  geometry  model  is  shown  in  Figure  19a.  The  geometry  included  a 
section  of  the  backward  facing  step  and  extended  approximately  three  step  heights  in  the  upstream 
and  downstream  directions  and  to  three  step  heights  above  the  slot  to  ensure  that  the  boundary 
conditions  do  not  influence  the  exhaust  velocity  distribution.  In  Figure  19a,  the  red  circles  indicate 
the  speaker  exits.  When  the  speakers  are  driven  by  a  single  frequency,  periodic  blowing  and  suction 
results  at  the  slot  exit.  The  peak  blowing  stroke  is  shown  in  Figure  19b,  where  the  color  represents 
the  wall  pressure  in  the  forcing  duct  and  the  velocity  is  shown  by  arrows  at  the  forcing  slot  exit.  The 
results  indicate  a  slight  spanwise  pressure  variation  (A p/ 1%)  from  the  center  of  the  chamber 
to  its  spanwise  edges.  The  exit  velocity  is  shown  to  be  essentially  uniform  (the  variations  seen  in 
the  figure  are  due  to  vectors  in  the  slot  boundary  layer). 

5.2.2.4  Open  loop  forcing  Numerous  open-loop  forced  simulations  were  performed  to  provide 
data  for  the  development  of  reduced  order  models  for  feedback  flow  control  (see  Section  5.2.3). 
These  open-loop  data  have  to  span  the  range  in  the  amplitude-frequency  parameter  space  that 
will  be  utilized  by  the  controller.  To  provide  these  data  sets,  simulations  with  different  forcing 
parameters  have  been  performed  and  analyzed. 

From  the  unforced  data  it  was  determined  that  the  vortical  structures  in  the  shear  layers  nat¬ 
urally  occur  at  Fn  ~  400Hz  at  x/H  =  2.  This  frequency  formed  the  basis  for  a  study  where  the 
blowing  and  suction  actuation  was  used  in  a  frequency  range  from  Ff  =  400Hz  to  Ff  =  1000Hz 
and  an  amplitude  ranging  between  A/ Uoo  =  0.01  and  A/Uoo  =  0.3,  resulting  in  the  time  dependent 
blowing  and  suction  velocity  u/(t )  =  A  si  n  (2tt  /Y  ) .  A  table  of  all  the  computed  cases  is  given  in 
Table  2;  representative  results  from  this  part  of  the  investigation  are  shown  in  this  section. 

When  forcing  is  applied  at  Ff  =  400Hz,  A/Uoo  =  0.1  (Figure  20),  the  density  spectra  taken 
at  the  five  downstream  locations  x/H  =  0,0.6, 1.3, 2, 2. 6  show  that  the  flow  initially  amplifies  the 
forcing  frequency  through  x/H  =  2.  At  x/H  =  2.6,  the  amplitude  starts  to  decay.  In  addition, 
the  first  harmonic  at  F  =  800Hz  is  amplified  between  x/H  =  1.3  and  x/H  =  2.  No  subharmonic 
frequency  is  discernable  in  the  data.  When  forcing  at  Ff  =  600Hz,  A/Uoo  =  0.1  (Figure  21),  the 
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U  velocity 


V  velocity 


Density 


(e)  (f) 

Figure  18:  POD  modes  and  time  coefficients  for  a),b)  u-velocity,  c),d)  v-velocity,  e),f)  density. 
Unforced  case. 
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(a)  Geometry  model  (b)  CFD  results 

Figure  19:  Blowing  and  suction  slot  analysis.  Pressure  distribution  in  the  blowing  suction  duct 
(color)  and  velocity  vectors  at  the  slot  exit. 


Frequency  (Hz) 


Figure  20:  Density  spectra  between  x/H  =  0  and  x/H  =  2.6,  F  =  400Hz,  A/Um  =  0.1. 


fundamental  is  strongly  amplified  between  x/H  =  0  and  x/H  =  0.6  and  decays  downstream.  In 
addition,  the  first  subharmonic  frequency  Fs  =  300Hz  is  amplified  as  well  and  decays  very  slowly 
downstream  of  x/H  =  1.3.  Finally,  forcing  at  Ff  =  800Hz,  A/Uoo  =  0.1  (Figure  22)  shows  the 
largest  amplitude  response  close  to  the  slot  of  the  three  forcing  cases  and  a  rapid  growth  of  the 
fundamental.  Downstream  of  x/H  =  0.6  it  decays  and  the  subharmonic  starts  to  develop,  indicating 
vortex  pairing  in  this  region. 

Performing  POD  on  the  density  results  in  the  modes  shown  in  Figure  23  for  a  forcing  frequency 
of  ff  =  400Hz.  In  the  figure,  POD  mode  isosurfaces  p  =  0.005kg/m  are  shown.  The  dominant 
modes  1  and  2  show  the  developing  shear  layer  vortex  street.  Modes  3  and  4,  which  develop  further 
downstream,  are  representative  of  the  spanwise  distortion  of  the  main  shear  layer  structures.  The 
time  coefficients,  Figure  23c,  clearly  show  the  forcing  frequency  in  modes  1  and  2. 

For  forcing  at  Ff  =  600Hz,  the  POD  modes  of  density  are  shown  in  Figure  24  (the  same 
isosurface  level  as  in  Figure  23  is  plotted).  The  development  of  the  POD  modes  is  similar  to 
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Figure  21:  Density  spectra  between  x/H  =  0  and  x/H  =  2.6,  F  =  600Hz,  A/Uoo  =  0.1. 


Figure  22:  Density  spectra  between  x/H  =  0  and  x/H  =  2.6,  F  =  800Hz,  A/Uoo  =  0.1. 
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(a)  Modes  1  and  2  (b)  Modes  3  and  4  (c)  POD  time  coefficients 

Figure  23:  POD  density  modes  and  time  coefficients,  F  =  400Hz,  AfUm  =  0.1. 


(a)  Modes  1  and  2  (b)  Modes  3  and  4  (c)  POD  time  coefficients 

Figure  24:  POD  density  modes  and  time  coefficients,  F  =  600Hz,  A/Uoo  =  0.1. 


the  Ff  =  400Hz  case,  although  the  emergence  of  the  modes  was  shifted  upstream  due  to  the  higher 
frequency,  in  accordance  with  theory  Ho  and  Huerre  [1984],  The  time  coefficients  exhibit  the 
same  behavior  as  in  the  previous  case,  but  the  subharmonic  character  of  modes  3  and  4  is  more 
pronounced.  Finally,  for  Ff  =  800Hz,  the  POD  results  are  shown  in  Figure  25.  The  figure  shows 
that  the  vortical  structures  develop  a  very  short  distance  downstream  of  the  step  (located  at  the 
inflow  boundary  of  the  box  shown  in  the  figures),  but  the  growth  saturates  quickly  and  the  highly 
coherent  structures  begin  to  show  three-dimensional  distortions  (as  indicated  by  POD  modes  3  and 
4)  further  downstream.  The  time  coefficients  show  that  modes  1  and  2  are  the  most  periodic  of  all 
investigated  forcing  cases  and  that  modes  3  and  4  oscillate  at  the  first  subharmonic  frequency.  In 
all  cases,  as  pointed  out  for  the  unforced  case,  the  density  POD  modes  represent  only  the  structures 
in  the  shear  layer  and  not  the  smaller  scale  motion  in  the  recirculation  region,  which  is  beneficial 
for  the  development  of  a  reduced  order  model  of  this  flow. 

The  optical  properties  of  the  flow  field  were  analyzed  using  the  beam  grid  described  above.  As 
for  the  unforced  case,  the  density  field  on  this  grid  was  integrated  from  the  wall  through  the  shear 
layer  and  the  effect  of  open  loop,  periodic  forcing  was  assessed.  Figure  26  shows  instantaneous 
plots  of  the  flow  structures  and  the  OPD.  There  is  a  strong  correlation  between  the  structures  in  the 
shear  layer  and  the  OPD  results,  as  was  observed  for  the  unforced  flow.  However,  when  forcing 
is  introduced,  the  structures  in  the  flow  exhibit  increased  spanwise  coherence  due  to  the  spanwise 
uniform  forcing. 

The  three-dimensional  simulations  described  above  all  showed  that  the  dominant  dynamics 
in  the  shear  layer  behind  the  backward  facing  step  are  essentially  two-dimensional,  at  least  in 
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(a)  Modes  1  and  2  (b)  Modes  3  and  4  (c)  POD  time  coefficients 


Figure  25:  POD  density  modes  and  time  coefficients,  F  =  800Hz,  A/U^  =  0.1. 


Figure  26:  Instantaneous  isosurface  of  flow  structures  and  associated  OPD  for  various  forcing 
conditions. 
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the  domain  investigated  in  this  research  project.  Especially  the  POD  analyses  showed  this  two- 
dimensionality  in  the  first  pair  of  modes,  which  are  representative  of  the  most  dominant  structures 
in  the  flow,  the  Kelvin-Helmholtz  vortices.  While  three-dimensional  effects  cannot  be  discarded 
for  a  complete  description  or  reconstruction  of  the  flow  field  from  the  POD  data,  it  is  posited 
that  for  the  development  of  feedback  flow  control  strategies,  a  two-dimensional  representation  of 
the  flow  captures  the  relevant  physical  processes.  This  allowed  for  drastically  reduced  simulation 
times  and  more  efficient  use  of  the  computational  time  available  to  the  project. 

To  develop  the  POD  mode  database  for  reduced  order  modeling  and  controller  design,  there¬ 
fore,  the  simulations  above  were  repeated  for  the  two-dimensional  grid.  The  results  obtained  from 
these  simulations  were  comparable  to  the  results  described  above.  Most  importantly,  the  natural 
shedding  frequency  matched  the  three-dimensional  results.  This  is  a  good  indicator  that  the  insta¬ 
bility  mechanism  that  results  in  the  Kelvin-Helmholtz  vortices  was  not  negatively  affected  by  the 
reduction  to  two  dimensions.  A  comparison  for  the  forced  simulations  showed  a  slightly  higher 
amplitude  of  the  fluctuations  for  the  two-dimensional  simulations.  This  was  expected  since  these 
simulations  could  be  viewed  as  ideal  in  the  sense  that  they  provide  perfect  spanwise  coherence. 
Another  way  to  think  of  these  results  is  that  the  structures  are  infinitely  long  in  the  spanwise  di¬ 
rection.  The  increased  coherence  of  the  structures  did  not  affect  the  initial  development  of  the 
shear  layer  targeted  with  feedback  control.  However,  further  downstream,  where  the  structures 
start  to  develop  spanwise  distortions  in  the  three-dimensional  simulations,  the  two-dimensional 
simulations  overpredicted  the  strength  of  the  structures. 

The  two-dimensional  simulations  were  performed  on  the  grid  shown  in  Figure  13.  Simulation 
results  for  the  unforced  or  open-loop  cases  will  not  be  shown  here;  results  for  the  validation  of  the 
feedback  control  strategy  are  presented  in  Section  5. 3.2. 2. 


5.2.3  Reduced  Order  Modeling 


5.2.3.1  Numerical  Reduction  The  procedure  in  section  4.2.2. 1  was  precisely  followed  and 
the  following  results  were  detected. 

The  quantity  of  greatest  concern  in  this  project  is  the  OPD,  which,  as  shown  in  Equations  57- 
59,  is  linearly  dependent  on  the  fluid  density  p(x,t).  Substituting  the  POD  decomposition  of  p 
into  Equation  57  yields 

m' 

n  =  1  +  Kgd  <Pi (x)  (6°) 

i=  1 

and  the  OPL  (Equation  58)  can  be  calculated  as 


OPL  =  [  n(s)ds  =  [ 

JL  JL 


1  +KGD'£ai(t)(pi(x) 


i=  1 


ds 


(61) 


OPL  =  L  +  KgdY  ai(t)  /  (pi(x)ds.  (62) 

Jl 

The  OPD  is  then  computed  from  Equation  59  such  that 


OPD  = 


m  r>  m  n 

L  +  KGDJ^ai(t)  Jl  (pi(x)ds  -  L  +  KGDJ^ai(t)  (pi(x)ds 


(63) 
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Simplifying  this  expression  gives 

m!  p  _  _ , 

OPD(x,t)  =  KGDY,ai(t)  /  (pi(x)  -  (pi(x)  ds  (64) 

t i  lJl  l  j  J 

where  the  OPD  spatial  modes  are  fL  ^<p,(x)  —  (Pi(x)'j  ds. 

Comparing  this  result  with  Equation  17  shows  that  the  OPD  POD  modes  can  directly  be  de¬ 
rived  from  the  density  POD  modes  by  using  the  definition  of  the  OPD.  This  is  important  because 
the  density  is  readily  available  from  the  CFD  simulations.  Furthermore,  because  of  this  direct  re¬ 
lationship,  the  minimization  of  the  time  coefficients  of  the  density  POD  modes  will  be  regarded  as 
the  main  control  goal  to  regulate  the  flow  to  reduce  OPD  fluctuations  using  feedback  flow  control. 


5.2.3.2  POD  Parameter  Study  A  large  parameter  study  was  carried  out  to  determine  an  ap¬ 
propriate  means  for  numeric  decomposition.  The  parameters  consisted  of  POD  vs.  DPOD,  spatial 
domain  size,  and  also  the  data  sets  to  be  used  in  the  decomposition.  The  spatial  domains  were 
defined  as 

XI  =  0<§<4  -1<£<1 

X2=  0<jj<2\  (65) 

X3—  l;j<)|-<2| 

The  three  spatial  domains  (XI,  X2,  X3)  adequately  contain  the  optical  aperture  of  interest,  whose 
center  is  located  at  x/H  =  2.  The  domain  size  study  helped  to  determine  if  limiting  the  amount  of 
information  in  the  POD  kernel  has  a  detrimental  effect  on  the  model. 

A  typical  forcing  input  for  a  given  frequency  and  amplitude  is  shown  in  Figure  27.  The  forcing 
begins  from  the  fully  developed  unforced  flow  computed  for  t  <  0.  The  first  five  cycles  of  forcing 
are  defined  as  the  opening  transient  in  which  the  flow  starts  to  react  to  the  forcing  signal.  The 
flow  then  locks  in  to  the  forcing  signal  for  the  remainder  of  the  duration  of  the  forcing.  When  the 
forcing  is  turned  off,  the  flow  undergoes  an  ending  transient  in  which  the  flow  shifts  back  into  its 
natural  state.  The  data  sets  used  for  the  POD/DPOD  study  were 

D1  =  Forcing 

D2  =  Forcing  +  U  n  forced  (66) 

D3  =  Forcing  +  U  n  forced  +  Starting  /  EndingT  ransients. 

The  resulting  POD/DPOD  model  was  then  validated  using  the  case  Ff  =  600Hz,  A/Uoo  =  0.05, 
which  had  been  removed  from  the  training  data  sets  for  model  validation  purposes.  To  quantify  the 
reconstruction  error,  the  root  mean  squared  error  of  the  local  and  instantaneous  density  field  was 
computed  as 

r - 1 1/2 

£rms(t)  =  [(p((Xi),t)-p((Xi),t))2  J  ,  (67) 

or,  expressed  as  a  percentage  of  the  mean  density. 


£(t) 


p((Xi)f) 


(68) 


The  parameters  for  the  dataset  study  are  shown  in  Table  3. 
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Figure  27:  Typical  forcing  input,  u(t\F,A),  for  an  open  loop  simulation.  Case  will  undergo  a  total 
of  four  defined  stages.  Starting  transient  (0  <  t  <  0.01),  Locked  in  (0.01  <t<  0.025),  Ending 
transient  (0.025  <t<  0.035)  and  unforced  (0.035  <  t  <  0.05). 


Table  3:  Summary  of  parameters  chosen  for  the  POD  dataset,  spatial  domain  and  method  study. 


Case 

Spatial  Domain 

Di 

Data  Set 
Ff  [Hz] 

A/Uoo  [%] 

e 

POD 

(0 

DPOD 

1 

XI 

D1 

400,  600,  800 

2.5,5,  10 

0.58% 

0.55% 

2 

X2 

Dl 

400,  600,  800 

2.5,5,  10 

0.67% 

0.52% 

3 

X3 

Dl 

400,  600,  800 

2.5,5,  10 

0.45% 

0.37% 

4 

XI 

D2 

400,  600,  800 

2.5,5,  10 

0.57% 

0.45% 

5 

X2 

D2 

400,  600,  800 

2.5,5,  10 

0.66% 

0.46% 

6 

X3 

D2 

400,  600,  800 

2.5,5,  10 

0.45% 

0.35% 

7 

XI 

D3 

400,  600,  800 

2.5,5,  10 

0.54% 

0.42% 

8 

X2 

D3 

400,  600,  800 

2.5,5,  10 

0.65% 

0.44% 

9 

X3 

D3 

400,  600,  800 

2.5,5,  10 

0.44% 

0.33% 
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The  POD  mode  sets  were  truncated  to  5  modes  for  POD  and  5x2  modes  for  DPOD.  Therefore, 
DPOD  was  expected  to  perform  better  at  reconstruction  of  the  density  field  than  POD.  The  error 
values  in  the  last  two  columns  in  Table  3  show  the  training  data  selection  and  spatial  domain  choice 
did  influence  the  reconstruction  error.  It  decreased  slightly  for  larger,  more  comprehensive  data 
sets.  Also,  the  reconstruction  error  decreased  proportionally  to  the  spatial  domain,  suggesting  that 
by  limiting  the  spatial  domain  size,  certain  flow  features  that  do  not  contribute  to  the  shear  layer 
physics  under  investigation,  such  as  the  recirculation  zone,  were  neglected.  Most  importantly,  the 
errors  indicate  that  it  is  crucial  to  retain  as  much  information  about  the  flow  field  as  possible,  shown 
by  the  fact  that  the  smallest  errors  were  obtained  using  dataset  D3,  which  includes  the  unforced 
and  forced  data  as  well  as  both  startup  and  shutdown  transients.  However,  the  main  outcome  of 
this  parameter  study  was  that  the  reconstruction  error  was  not  a  good  way  to  quantify  which  mode 
set  to  choose  as  the  final  numeric  model.  All  of  the  error  values  were  acceptable  (<  1  per  cent), 
which  would  suggest  that  all  these  parameter  combinations  would  be  adequate. 

Scrutinizing  the  differences  between  POD  and  DPOD  provides  insight  into  the  shear  layer 
dynamics  as  forcing  is  applied  to  the  flow.  From  theory  as  well  as  the  experimental  data  obtained 
in  this  project,  it  is  well  known  that  the  shear  layer  is  extremely  susceptible  to  periodic  forcing. 
Due  to  the  flow’s  instability,  small  perturbations  over  a  wide  frequency  range  are  amplified  and 
result  in  Kelvin-Helmholtz  vortices  [see  e.g.  Oster  and  Wygnanski,  1982],  The  shift  modes  for 
the  DPOD  mode  sets  ((pi, 2),  shown  in  Figures  28,  31,  and  34,  model  the  transient  change  from 
the  natural  shedding  to  a  forced  state,  but  the  shift  modes  show  a  change  in  wavelength  compared 
their  corresponding  main  mode,  indicating  that  there  was  not  a  slow  shift  from  the  natural  to  a 
forced  state.  The  data  suggested  that  for  the  shear  layer,  the  flow  response  was  different  because 
the  frequency  band  in  which  lock-in  occurs  is  much  larger  than  for  flows  such  as  the  cylinder  wake. 
As  seen  from  these  DPOD  modes,  the  shear  layer  structures  assume  a  different  wave  length  when 
forced  at  a  given  frequency.  The  transient  behavior  is  extremely  fast,  thus  making  the  underlying 
concept  of  DPOD  questionable.  In  addition,  the  spatial  modes  for  the  DPOD  decomposition  lack 
physical  relevance. 

The  POD  models  for  cases  1,  6,  and  8  are  shown  in  Figures  29,  32  and  35,  respectively.  The 
mode  sets  show  the  mean  flow  and  the  first  2  mode  pairs.  As  expected,  there  is  a  distinct  size/wave¬ 
length  change  for  the  2  mode  pairs.  The  modes  in  the  POD  mode  sets  also  look  physically  viable, 
unlike  the  DPOD  mode  sets.  Therefore,  POD,  not  DPOD,  was  chosen  as  the  preferred  modeling 
approach  for  the  shear  layer. 

When  the  reconstruction  error  is  evaluated  as  a  function  of  time,  as  in  Equation  67,  DPOD 
should  better  represent  the  flow  field  in  the  transient  regions  whereas  POD  should  fail  to  model 
these  transitions  from  one  flow  state  to  another.  Figures  30,  33,  36  show  that  this  was  not  the  case. 
In  fact,  the  reconstruction  error  is  smaller  for  the  POD  approximation  than  for  the  DPOD  one.  This 
corroborates  the  assertion  that  the  five  supplemental  shift  modes  were  not  really  modeling  transient 
flow  effects  in  the  shear  layer. 
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Figure  28:  Case  1:  DPOD  modes  <p,/  for  parameters  shown  in  Table  3.  Left  column:  Main  modes 
(pn,  right  column:  Shift  modes  <Pq. 
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Figure  29:  Case  1:  POD  modes  <p,  for  parameters  shown  in  Table  3.  a)  Mean  flow  mode,  b)  and  c) 
first  fluctuating  mode  pair,  d)-e)  first  harmonic  mode  pair. 
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Time  [s] 

Figure  30:  Case  1 :  (Top)  Reconstruction  error  as  a  function  of  time  of  numerical  model,  (-)  POD, 
( — )  DPOD.  (Bottom)  Forcing  signal  for  validation  case.  Transient  period  in  flow  field  begins 
around,  lock  in  region  to  t  =  0.03 s,  ending  transient  until  t  =  0.03 5s  where  the  natural  flow  state 
occurs. 
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Figure  31:  Case  6:  DPOD  modes  (pij  for  parameters  shown  in  Table  3.  Left  column:  Main  modes 
(pn,  right  column:  Shift  modes  (pa. 
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Figure  32:  Case  6:  POD  modes  (pi  for  parameters  shown  in  Table  3. 

The  only  parameter  left  to  choose  was  the  spatial  domain.  The  spatial  domain  sets  were  chosen 
to  be  the  entire  flow  field  behind  the  step  (XI),  just  the  shear  layer  neglecting  the  recirculation 
zone  (X2)  and  the  shear  layer  over  the  optical  aperture  (X3).  Because  the  reconstruction  error,  as 
shown  in  Table  3,  did  not  provide  a  reliable  criterion  for  which  spatial  domain  was  appropriate, 
the  domain  was  chosen  on  a  physical  basis.  The  first  spatial  domain,  XI,  limited  the  ability  to 
capture  the  dynamics  of  interest  by  retaining  undesirable  flow  physics,  such  as  the  recirculation 
zone,  in  the  domain.  As  shown  in  Figure  37,  the  mode  amplitudes  obtained  for  the  validation  case 
( Ff  =  600Hz,  A/Uoo  =  0.10)  showed  some  low  frequency  content,  which  originated  most  likely 
from  the  recirculation  zone.  Figure  38,  the  corresponding  data  for  domain  X2,  exhibits  much  more 
periodic  mode  amplitudes,  which  are  representative  of  the  structures  in  the  shear  layer. 
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Time  [s] 


Figure  33:  Case  6:  (Top)  Reconstruction  error  as  a  function  of  time  of  numerical  model,  (-)  POD, 
( — )  DPOD.  (Bottom)  Forcing  signal  for  validation  case.  Transient  period  in  flow  field  begins 
around  t  =  0.005,?,  lock  in  region  to  t  =  0.03?,  ending  transient  until  t  =  0.035?  where  the  natural 
flow  state  occurs. 
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Figure  34:  Case  8:DPOD  modes  (<p,y)  for  parameters  shown  in  Table  3.  (a)(c)(e)(g)(i)  Main 
inodes(<p/i,  (b)(d)(f)(h)(k)  Shift  modes  (<p;- 2). 
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Figure  35:  Case  8:  POD  modes  (<p;)  for  parameters  shown  in  Table  3. 
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Figure  36:  Case  8:  (Top)  Reconstruction  error  as  a  function  of  time  of  numerical  model,  (-)  POD, 
( — )  DPOD.  (Bottom)  Forcing  signal  for  validation  case.  Transient  period  in  flow  field  begins 
around  t=0.005s,  lock  in  region  to  t=  0.003s,  ending  transient  until  t=0.0035s  where  the  natural 
flow  state  occurs. 
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Figure  37:  Case  1:  POD  mode  amplitudes  (a,-)  for  parameters  shown  in  Table  3. 
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(a)  Mode  a\ 
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Figure  38:  Case  8:  POD  mode  amplitudes  (a,-)  for  parameters  shown  in  Table  3. 
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5.3  System  Identification 

Another  important  parameter  was  the  size  of  the  training  data  set  for  the  model.  The  number 
and  span  of  training  cases  for  the  WNARX  model  is  presented  in  Table  4.  A  total  of  12  open 
loop  cases,  all  of  which  contained  starting  and  ending  transients  from  the  unforced  flow  state  and 
back  to  it,  were  computed  to  understand  the  influence  of  actuation  with  varying  frequency  and 
amplitude  on  the  flow.  The  results  of  this  investigation  showed  that  the  time  coefficients  reacted 
almost  linearly  to  the  blowing  and  suction  amplitude,  i.e.  the  response  of  the  mode  amplitudes, 
aj(t),  scaled  linearly  with  amplitude  input.  In  contrast,  the  flow  response  was  highly  nonlinear 
with  respect  to  the  forcing  frequency.  Thus,  the  three  training  data  sets  highlighted  in  Table  4  were 
chosen  to  provide  a  basis  space  for  the  WNARX  model.  The  case  Ff/Fn  =  l,  A/Uoo  =  0.10  was 
chosen  to  be  the  validation  case  for  the  model.  A  summary  of  final  parameters  for  the  dynamic 
model  is  presented  in  Table  5. 


Table  4:  Summary  of  cases.  / :  WN  training  cases,  o:  validation  case. 


400Hz 

600Hz 

800Hz 

A/Uoo  =  0.1 

/ 

0 

/ 

A/Uoo  =  0.05 

X 

/ 

X 

A/t/oc  =  0.025 

X 

X 

X 

A/Uoo  =  0.0125 

X 

X 

X 

The  WNARX  model  was  validated  for  an  off  design  flow  case  for  which  the  forcing  signal  was 
turned  on  at  simulation  time  t  =  0s,  at  which  point  the  flow  goes  through  a  transient  period  before 
locking  into  the  forcing  frequency.  The  forcing  was  then  turned  off  at  t  =  0.025s  (corresponding  to 
15  forcing  periods)  to  reestablish  the  unforced  flow  state  (see  Figure  27).  As  seen  in  Figure  39,  the 
model  captures  the  lock-in  region  of  the  periodic  forcing  very  well.  Once  the  forcing  was  turned 
off  at  t  =  0.025s,  the  model  accurately  predicted  the  type  of  nonlinear  signal  in  the  unforced  flow. 
Expecting  an  exact  replication  of  the  unforced  time  coefficients  is  unrealistic  since  the  signal  was 
extremely  nonperiodic.  However,  the  important  point  is  that  the  model  of  the  unforced  flow  does 
not  decay  to  zero  over  time.  This  indicates  that  there  is  a  periodic  attractor  to  the  nonlinear  function 
for  the  WNARX  system.  The  similarities  in  periodic  trends  furthermore  suggest  that  the  attractor 
is  near  the  solution  of  the  unforced  state. 

With  the  development  of  this  model,  the  feedback  control  problem  of  the  shear  layer  behind  a 
backward  facing  step  had  been  transformed  into  the  problem  of  designing  a  controller  for  the  POD 
time  coefficients.  At  this  point,  the  previous  definition  of  the  control  goal,  namely  minimizing  the 

Table  5:  Summary  of  parameters  chosen  for  the  WNARX  model. 

Mode  at  Wavelets  Regressors 
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Figure  39:  Off  design  validation  of  the  four  mode  WNARX  model  for  flow  case  of  Ff  =  600Hz 
and  A/Uoo  =  0. 1 .  WNARX  output  (-),  POD  model  (-  -). 
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Table  6:  Summary  of  parameters  chosen  for  surface  mounted  sensor  estimation  scheme. 


Study 

#  of  Sensors 

A(£) 

Data  Set 

Pressure  SPOD  TC 

LSE 

Error  e(t ) 
ANNE  WNE 

1-3 

47 

2.5%  5%  10% 

/ 

- 

21% 

19% 

12% 

4-6 

47 

2.5%  5%  10% 

- 

/ 

22% 

20% 

12% 

7-9 

47 

2.5%  5%  10% 

/ 

/ 

11% 

9.5% 

4.6% 

10-12 

47 

5%  10% 

/ 

/ 

10.9% 

9.6% 

4.6% 

13-15 

8 

5%  10% 

/ 

/ 

11.2% 

10.8% 

4.7% 

16-18 

6 

5%  10% 

/ 

/ 

18.8% 

16.5% 

5.2% 

19-21 

4 

5%  10% 

/ 

/ 

20.5% 

19.2% 

8.1% 

22-24 

2 

5%  10% 

/ 

/ 

23.1% 

19.9% 

21% 

density  fluctuations  in  the  flow  field,  was  replaced  by  a  much  more  tractable  problem:  Design  a 
controller  for  the  model  of  the  POD  time  coefficients  whose  dimensionality  is  orders  of  magnitude 
less  than  that  of  the  underlying  CFD  flow  field  data.  Finding  a  model  that  described  a  flow  field  that 
was  not  included  in  the  model  design  with  the  fidelity  shown  was  a  major  step  toward  developing 
successful  feedback  flow  control  strategies  for  the  free  shear  layer  flow. 

5.3.1  Feedback  Control 

5.3.1. 1  State  Estimation  A  parameter  study  was  conducted  for  the  state  estimator  formulation. 
Because  the  density  of  the  fluid  cannot  be  directly  measured  on  the  surface  behind  the  step,  pressure 
variations,  which  directly  correlate  to  the  density  fluctuations  in  the  flow  (with  the  assumption  of 
constant  temperature),  were  chosen  as  viable  surface  measurements.  Five  forcing  cases,  which 
contained  starting  and  ending  transients  (see  Table  2)  were  used  as  state  estimator  training  data. 
The  case  Ff  =  600Hz,A/t/oo  =  0.1  (Table  2)  was  reserved  for  validation  purposes.  It  should  be 
noted  that  the  error  of  the  training  data  is  bounded  above  by  the  error  of  the  validation  case,  that  is 

\\^T rain  \ \  f  ||^Va/||j  (69) 

for  all  cases.  Therefore,  the  error  of  the  validation  case  will  suffice  to  determine  the  performance 
of  the  estimation  method.  Table  2  shows  that  the  estimator  was  interpolating  between  cases;  the 
accuracy  of  the  estimator  outside  the  training  region  has  not  been  verified. 

Determining  the  appropriate  parameters  for  the  state  estimator  began  by  defining  the  physical 
locations  of  the  sensors.  The  floor  behind  the  step  from  0<xs  <  2.5 H  contained  a  total  of  47  pos¬ 
sible  sensor  locations.  All  three  estimation  methods  (LSE,  ANNE,  WNE)  were  applied  to  the  full 
state  sensor  array  (i.e.  dim{xs )  =  47)  to  determine  the  best  performance  of  the  estimation  methods. 
The  sensor  array  was  then  down-sampled  to  the  minimum  number  of  sensors  (i.e.  dim(xs )  =  2)  and 
incrementally  increased  until  the  error  converged  to  the  full  state  sensor  estimation  performance. 
The  parameters  of  the  training  methods  for  ANNE  and  WNE  were  held  constant  throughout  this 
study.  A  total  time  history  of  25  time  steps  was  used  in  the  formulation  of  the  regression  vector 
(i.e.  n  =  25  in  Equation  38).  The  results  are  summarized  in  Table  6. 

Figure  40a)  shows  the  error  in  the  estimation  of  mode  1  as  a  function  of  the  number  of  sensors. 
The  results  indicate  that  all  of  the  methods  rapidly  converge  to  approximately  their  final  perfor¬ 
mance  level  with  the  use  of  8  sensors  spaced  equally  between  x/H  =  0.25  and  x/H  =  2.  More 
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Figure  40:  Estimation  error  for  validation  case  Fj  =  600Hz,  A/Uoa  =  0.1.  a)  Root  mean  squared 
error  of  mode  1,  a\ (t),  all  estimators,  b)  Root  mean  squared  error  of  modes  1-4  using  the  wavenet 
estimator  (WNE). 

interestingly,  the  wavenet  estimation  method  resulted  in  only  half  the  error  of  the  other  methods 
( 4.5  per  cent  for  dim(xs )  =  8).  In  Figure  40b),  the  errors  for  all  four  modes  as  computed  using  the 
WNE  estimator  are  plotted  as  a  function  of  the  number  of  sensors.  The  plot  indicates  that  while  the 
error  increased  somewhat  for  the  higher  modes,  all  modes  were  converged  when  using  only  eight 
sensors. 

The  results  obtained  with  the  full  state  sensor  array  along  with  surface  POD  time  coefficients 
(Study  10-12)  are  shown  in  Figure  41.  This  is  the  best  possible  result  given  the  amount  of  infor¬ 
mation  and  infinite  training  time.  The  conclusion  of  this  study  represented  a  trade  off  between  the 
number  of  sensors  needed  for  precise  estimation  and  the  wavenet  training  time.  The  goal  was  to 
have  the  minimal  number  of  sensors  for  accurate  estimation  while  maintaining  a  physically  feasible 
sensor  configurations. 

Figure  42  shows  the  estimation  results  for  only  two  sensors  in  the  sensor  configuration.  From 
the  error  computations,  it  is  clear  that  the  estimators  need  larger  sensor  arrays  for  accurate  estima¬ 
tions  of  the  flow  field.  Figure  43  is  the  optimal  sensor  configuration  which  was  determined  to  be  an 
array  of  eight  sensors  between  xs/H  =  0.25  and  x/H  =  2.  RMS  errors  were  on  the  order  of  5  per 
cent  for  this  sensor  configuration,  which  was  equivalent  to  the  error  of  the  estimation  using  the  full 
state  sensor  estimate.  Figure  44  shows  a  comparison  of  the  actual  time  coefficients  computed  from 
Equation  17  with  the  simulated  WNE  computed  from  Equation  27  using  the  eight  sensors.  The 
estimator  captures  both  the  phase,  frequency,  and  amplitude  of  the  flow  states  for  the  validation 
case.  At  this  point  the  density  field  could  be  reconstructed  with  an  error  of  less  than  5  per  cent  of 
the  original  flow  field  using  only  eight  sensors  by  combining  surface  POD  (Equation  17)  and  the 
flow  state  estimate  (Equation  27)  within  the  forcing  parameter  space. 

5.3.2  Adaptive  Control 

Direct  adaptive  feedback  control  [see  Fagley  et  al.,  2009,  for  more  details]  was  chosen  to  close  the 
feedback  loop.  Adaptability  allows  for  uncertainties  when  scaling  the  controller  for  validation  in 
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Figure  4 1 :  Estimation  using  full  sensor  array  of  forcing  validation  case  (Ff  =  600Hz,  A/U0o  =  0.1). 
Forcing  on  for  first  half  of  simulation,  off  for  second,  a)  Study  10:  LSE  estimations  (e  =  10.9%) 
b)  Study  11:  ANNE  estimations  (e  =  9.6%)  c)  Study  12:  WNE  estimations  (e  =  21%). 
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Figure  42:  Estimation  using  smallest  sensor  array  of  forcing  validation  case  (Ff  =  600Hz,  A/Uoo  = 
0.1).  Forcing  on  for  first  half  of  simulation,  off  for  second,  a)  Study  22:  LSE  estimations  (£  = 
23.1%)  b)  Study  23:  ANNE  estimations  (£  =  19.9%)  c)Study  24:  WNE  estimations  (£  =  4.6%). 
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Figure  43:  Estimation  using  full  sensor  array  of  forcing  validation  case  (Ff  =  600Hz,  A/Uoo  =  0.1). 
Forcing  on  for  first  half  of  simulation,  off  for  second,  a)  Study  13:  LSE  estimations  (£  =  11.2%) 
b)  Study  14:  ANNE  estimations  (£  =  10.8%)  c)  Study  15:  WNE  estimations  (£  =  4.6%). 
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Figure  44:  Comparison  of  first  four  POD  time  coefficients  for  the  validation  case  Ff  =  600Hz, 
A/Uoo  =  0.1.  — ,  CFD  results, - ,  WNE  with  eight  sensors.  Forcing  was  active  for  t  <  0.025s. 


CFD  simulations  and  experiments.  The  basic  equations  describing  direct  adaptive  control  are 


U  =  Geey 
Ge  =  —eyey  ye , 


(70) 


where  Ge  is  the  gain  matrix,  ye  is  the  adaptability  weight,  and  ey  is  the  error  between  output  and 
desired  reference  signal, 

ey  =  a  —  aref.  (71) 

For  multi  input  multi  output  (MIMO)  systems,  ey  and  ye  are  matrices  of  size  nout  x  n,M.  Also,  the 
gain  matrix  is  of  size  n,w  x  nout.  The  derivative  must  be  approximated  numerically,  because  no 
analytic  solution  exists.  Here,  the  fourth  order  Adams-Bashforth  method, 


Ten+ 1 


-  G, 


en 


At 


=  e 


'V  H - Vev  .  -| - 1 - £y  o 

yn  2  sn—l  yn—2  g  yn—S 


5 

12 


(72) 


was  utilized  to  determine  the  gain  matrix  derivative.  The  feedback  parameters  associated  with 
this  control  strategy  are  primarily  which  mode  is  used  for  feedback  and  the  adaptability  weights 
which  are  typically  less  than  one.  Stability  of  this  type  of  control  system  is  only  proven  for  linear 
systems  Fuentes  and  Balas  [2000].  Stability  margins  cannot  be  shown  for  our  nonlinear  system  of 
equations  with  adaptive  control;  however,  stable  simulations  provide  empirical  evidence. 


5.3.2. 1  Feedback  control  of  the  WNARX  model  Developing  the  components  for  a  closed 
loop  simulation  is  a  multi-step  iterative  process.  The  model  developed  above  provides  accurate 
predictions  of  the  mode  amplitudes  when  the  flow  is  forced  within  a  given  frequency  and  parameter 
range,  including  starting  and  ending  transients.  However,  it  remains  to  be  seen  if  the  model  is 
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Figure  45:  Block  diagram  of  closed  loop  simulation  strategy.  Step  1:  Design  a  controller  for 
WNARX  ROM.  Step  2:  Apply  controller  to  CFD  simulation  with  state  estimator  to  verify  model 
adequately  predicts  closed  loop  dynamics. 


capable  of  adequately  simulating  the  highly  nonlinear  dynamics  expected  for  the  closed  loop  case. 
Furthermore,  a  feedback  scenario  allows  for  a  parameter  study  of  adaptive  control  algorithms. 
The  WNARX  model  allows  for  very  quick  simulation  times,  so  that  a  parameter  study  can  be 
carried  out  very  quickly.  The  parameters  were  adjusted  to  feed  back  different  combinations  of 
modes  and  their  derivatives  along  with  preconditioned  adaptability  weights.  Once  desirable  results 
were  achieved  with  the  model  in  a  closed  loop  simulation,  the  designed  control  algorithm  with  the 
corresponding  feedback  mode  combination  and  weights  was  scaled  up  and  applied  in  a  closed  loop 
CFD  simulation  to  validate  both  the  WNARX  ROM  system  and  the  adaptive  controller.  A  diagram 
of  the  two  parts  of  this  approach  is  shown  in  Figure  45. 

After  the  parameter  study  discussed  above,  it  was  determined  that  the  POD  mode  amplitude  a\ 
and  its  time  derivative,  a\ ,  were  the  best  parameters  to  be  regulated  in  the  feedback  control  system. 
The  derivative  of  a\  was  computed  by  an  implicit  Euler  approximation.  Because  this  can  be  a  poor 
approximation  of  the  derivative  and  its  susceptibility  to  noise,  a  moving  average  filter  was  added 
to  smooth  the  estimated  derivative.  The  initial  idea  was  that  the  OPD  would  be  reduced  by  simply 
reducing  the  mode  amplitudes.  Feedback  of  states  a\  and  ci\  allowed  for  excellent  controllability 
of  the  mode  amplitudes  as  shown  in  Figure  46.  Note  that  by  controlling  mode  1,  mode  2  was 
controlled  as  well  because  these  modes  represent  the  traveling  wave  character  of  the  shear  layer 
structures. 

In  this  simulation,  the  open  loop  forced  flow  was  used  as  the  initial  condition  for  the  closed  loop 
simulation  to  create  periodicity  in  the  flow  and  to  improve  startup  performance  of  the  controller 
when  the  loop  was  closed.  Figure  46  shows  the  time  coefficients  for  the  four  mode  model.  Periodic 
forcing  was  applied  for  t  <  0.015s,  at  which  point  the  closed  loop  control  was  switched  on  for  a 
time  period  of  0.015s  <  t  <  0.035s,  when  the  control  is  turned  off  and  unforced  flow  redeveloped 
for  t  >  0.035s.  As  shown  in  the  figure,  the  controller  performs  well,  reducing  the  amplitudes  of 
the  time  coefficients  to  approximately  35  per  cent  of  the  unforced  state.  As  a  final  step  to  verify 
the  efficacy  of  this  control  approach  for  aero-optical  problems,  the  density  field  was  reconstructed 
using  the  closed  loop  simulation  results  of  the  time  coefficients,  shown  in  Figure  46,  and  their 
corresponding  spatial  modes  (Figure  32).  The  reconstructed  density  field,  p(x,y,t),  allows  for  the 
evaluation  of  the  effect  of  the  three  different  forcing  scenarios  (unforced,  open-loop  forced,  and 
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Figure  46:  Feedback  results  using  adaptive  feedback  control.  Periodic  forcing  for  0s  <  t  <  0.0155, 
closed  loop  simulation  for  0.0155  <  t  <  0.0355,  and  unforced  simulation  for  t  >  0.0355. 
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Figure  47:  Feedback  results  using  adaptive  feedback  control.  Periodic  forcing  for  0s  <  t  <  0.015s, 
closed  loop  simulation  for0.015s<t<0.035s,  and  unforced  simulation  for  t  >  0.035s.  Verti¬ 
cal  solid  lines  indicate  contours  of  p'(y,t )  normalized  by  the  maximum  density  fluctuation.  The 
Horizontal  line  shows  the  maximum  density  fluctuation  at  a  given  time. 


feedback  controlled),  especially  their  effect  on  the  density  fluctuations.  Because  the  density  field 
is  three  dimensional  in  (x,y,t)-space,  it  is  difficult  to  visualize  the  flow  field  dynamics.  Here,  the 
standard  deviation  of  the  density  field  was  computed  for  visualization  using 


Nx 


Nxt 


- X 

P(xi,y,t)-p(x,y,t) 


(73) 


Note  that  the  mean  was  taken  in  the  x-direction.  In  effect,  with  these  two-dimensional  simulations, 
perfect  spanwise  coherence  was  assumed.  In  Figure  47,  contours  of  p'(y )  at  discrete  times  are 
plotted  as  vertical  lines.  The  figure  shows  that  the  magnitude  of  p'(y)  as  well  as  the  extent  of 
the  distortions  in  the  y-direction  were  significantly  reduced  when  feedback  control  was  active.  In 
addition,  max(p/(y))  is  plotted  for  all  times.  It  corroborates  the  reduction  of  the  density  fluctuations 
for  the  feedback  controlled  flow  field. 

As  a  final  performance  metric,  the  OPD  for  a  beam  passing  through  this  flow  field  was  com¬ 
puted  using  Equations  (57)-(59).  For  the  2D  simulations,  the  aperture  size  was  1.5  <  x/H  <  2.5 
with  unit  width.  The  OPD  at  the  point  of  interest,  x/H  =  2,  is  plotted  in  Figure  48,  which  shows 
that  the  OPD  was  drastically  reduced  during  the  closed  loop  portion  of  the  simulation,  both  in 
comparison  to  the  periodically  forced  flow  and  to  the  unforced  flow. 


5. 3. 2. 2  Feedback  control  in  the  CFD  simulation  The  final  validation  of  the  controller  devel¬ 
oped  during  this  research  effort  was  performed  by  implementing  the  controller  in  the  Cobalt  CFD 
simulations.  Hooks  had  been  added  in  the  Cobalt  CFD  code  during  an  earlier  AFOSR  funded 
STTR  project  between  Cobalt  Solutions,  LLC,  and  the  US  Air  Force  Academy.  These  hooks 
make  sensor  information  available  to  Matlab®  which  handles  the  controller  computations.  Af¬ 
ter  the  actuator  output  had  been  determined,  it  was  passed  back  to  the  Cobalt  simulation  using 
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Figure  48:  Calculated  optical  path  difference  at  x/H  =  2,  y/H  =  0  for  the  reconstructed  flow  field 
of  the  closed  loop  simulation.  Periodic  forcing  for  0s  <  t  <  0.015s,  closed  loop  simulation  for 
0.015s  <t<  0.035s,  and  unforced  simulation  for  t  >  0.035s. 

externally  controlled  blowing  and  suction  boundary  conditions.  The  program  communication  be¬ 
tween  Matlab®  and  Cobalt  allows  for  great  flexibility  when  incorporating  feedback  control,  filters, 
state  estimators,  etc.  within  a  CFD  simulation.  More  importantly,  it  allows  for  utilization  of  the 
exact  programs  developed  in  the  previous  sections,  which  significantly  reduces  the  possibility  of 
program  errors. 

For  the  first  validation  step,  the  controller  in  the  previous  section  was  directly  used  in  the  CFD 
simulation  in  conjunction  with  the  state  estimator  developed  in  Section  5. 3. 1.1.  The  feedback 
controlled  simulation  proceeded  as  follows:  First,  the  Cobalt  simulation  was  advanced  one  time 
step.  The  new  data  at  the  sensor  locations  (predetermined,  see  Section  4.2.4.2)  was  then  passed 
to  the  Matlab®  state  estimator  to  estimate  the  POD  mode  amplitudes;  the  estimation  was  seen  to 
be  essentially  the  same  as  what  the  model  predicted.  These  POD  mode  amplitude  estimates  were 
then  input  into  the  control  algorithm,  whose  output  was  converted  to  a  blowing  and  suction  mass 
flow  rate  for  the  blowing  and  suction  slot.  Finally,  this  information  was  passed  back  to  Cobalt  as  a 
new  boundary  condition  value  to  be  used  in  the  subsequent  CFD  iteration. 

After  completing  the  simulation  with  this  controller,  the  density  field  data  was  used  to  compute 
the  OPD.  First,  as  before,  the  density  fluctuations  were  computed  and  plotted  using  the  same 
method  as  for  the  WNARX  validation  (Figure  47).  Figure  49  shows  that  the  controller  developed 
using  the  ROM  had  a  pronounced  effect  in  reducing  the  density  fluctuations,  similar  to  the  effect 
observed  in  the  model  simulation  results. 

The  OPD  results  in  Figure  50  show  that  the  controller  (active  for  t  >  0.025s)  reduces  the  OPD, 
but  the  reduction  was  slower  than  predicted  by  the  WNARX  model.  This  was  most  likely  due  to 
discrepancies  between  the  reduced  order  model  and  the  CFD  simulation  results,  indicating  that 
the  ROM  did  not  quite  capture  all  the  intricate  nonlinear  dynamics  of  the  flow  field  which  were 
resolved  in  the  CFD  simulation.  As  shown  in  Table  5,  the  model  assumed  that  modes  a\  and  a2  are 
completely  decoupled  from  modes  a 3  and  a 4.  This  was  likely  the  most  dramatic  modeling  shortfall 
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Figure  49:  Density  analysis  of  CFD  closed  loop  simulation.  Periodic  forcing  for  0s  <  t  <  0.025s, 
closed  loop  simulation  for  t  >  0.025s.  Vertical  solid  lines  indicate  contours  of  p'iy.t)  normalized 
by  the  maximum  density  fluctuation.  The  Horizontal  line  shows  the  maximum  density  fluctuation 
at  a  given  time. 


Figure  50:  OPD  as  a  function  of  time.  CFD  results  with  WNARX  controller  directly  substituted 
into  CFD  closed  loop  simulation. 
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and  the  reason  that  the  model  did  not  correctly  capture  the  nonlinear  dynamics  in  the  closed  loop 
simulation.  Further  analysis  of  this  closed  loop  simulation  was  performed  to  better  understand 
these  differences. 

Figure  50  shows  that  it  took  approximately  At  =  0.025s  after  the  controller  activation  for  the 
control  to  be  effective  in  reducing  the  OPD.  Full  control  was  achieved  for  t  >  0.05s,  at  which 
time  the  large  amplitude  oscillations  due  to  flow  transients,  have  convected  downstream.  The 
OPD  amplitude  was  reduced  by  approximately  50  per  cent.  Scrutinizing  Figure  51e,  the  time 
history  of  the  control  input,  corroborates  that  at  time  t  =  0.05s,  the  controller  had  started  to  achieve 
the  control  goal  of  minimizing  the  POD  mode  1  amplitude  (Figure  51a)  and  reduced  the  forcing 
amplitude  to  approximately  A/Uoo  —  0.01,  which  was  less  than  one  per  cent  of  the  free  stream 
velocity. 

The  remaining  POD  mode  time  coefficients  are  shown  in  Figure  51b-d.  Since  POD  mode  2 
(Figure  51c)  is  the  complement  to  mode  1  to  comprise  a  traveling  wave,  it  was  not  surprising  that 
its  amplitude  was  reduced  in  unison  with  mode  1 .  However,  modes  3  and  4  behaved  differently;  the 
main  effect  of  forcing  on  these  modes  seemed  to  be  a  stabilization  of  their  oscillation  frequencies 
and  also  their  amplitudes. 

When  this  research  project  was  started,  it  was  believed  that  a  reduction  of  the  OPD  would  have 
to  be  coupled  to  the  minimization  of  these  mode  amplitudes,  since  each  mode  pair  represents  they 
flow  state  created  when  the  shear  layer  is  forced  with  a  given  frequency.  However,  the  results 
indicate  that  the  shear  layer  is  far  too  unstable  and  quickly  moves  away  from  the  natural  periodic 
attractor  when  forcing  is  applied.  A  comparison  of  the  OPD  results  with  the  POD  mode  amplitudes 
seems  to  suggest  that  the  desired  flow  state  is  in  fact  achieved  by  introducing  a  new  periodic  state 
which  reduces  the  OPD  for  a  given  aperture  location  and  size.  Because  discrepancies  between  the 
reduced  order  model  and  the  CFD  simulations  did  exist,  in  the  final  step  of  this  research  project, 
the  controller  parameters  in  Equation  70  were  adjusted  by  scrutinizing  the  CFD  results  directly 
to  increase  closed  loop  performance,  efficacy  and  efficiency.  These  results  are  presented  in  the 
following. 

Since  it  was  determined  from  the  initial  feedback  controlled  CFD  results  that  information  about 
POD  modes  3  and  4  needed  to  capture  the  transient  behavior  better,  a  combination  of  mode  ampli¬ 
tudes  and  their  derivatives  were  fed  back  in  numerous  CFD  simulations  to  determine  the  optimum 
combination  for  the  adaptive  control  algorithm.  It  was  found  that  feeding  back  the  time  coefficients 
of  one  of  the  modes  of  the  next  mode  pair,  a 3,  with  an  aggressive  adaptability  weight,  ye  =  1,  intro¬ 
duced  this  periodic  attractor,  which  effectively  reduced  the  density  fluctuations  and  therefore  the 
OPD  over  a  given  aperture.  Figures  52  and  53  show  the  final  successful  closed  loop  simulation  and 
the  corresponding  reduction  of  the  OPD  to  approximately  30  per  cent  of  its  original  value.  This 
presents  a  performance  improvement  of  almost  50  per  cent  over  the  original  controller.  In  addition, 
the  fluctuation  amplitude  was  drastically  reduced  when  compared  to  the  periodically  forced  flow  as 
well  as  when  compared  to  the  unforced  flow.  Scrutinizing  the  control  signal,  it  was  observed  that 
the  controller  introduces  two  harmonic  frequencies  into  the  flow,  the  lower  of  which  was  approxi¬ 
mately  Ff  ~  720Hz.  Interestingly,  the  control  amplitude  did  not  decline  as  the  controller  became 
effective,  as  initially  anticipated.  In  contrast,  the  amplitude  seemed  to  stabilize  at  A/Uoo  —  0.04, 
which  was  larger  than  for  the  controller  obtained  directly  from  the  ROM. 

The  results  from  this  closed  loop  simulation  supported  the  idea  that  excitation  of  frequencies 
that  are  unstable  further  upstream  (closer  to  the  step)  has  a  much  larger  effect  on  the  OPD,  even 
if  the  aperture  is  located  downstream,  than  forcing  at  a  frequency  that  is  commensurate  with  the 
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Figure  5 1 :  Mode  amplitudes  and  control  output  of  CFD  closed  loop  simulation. 
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naturally  occurring  frequency  at  the  aperture  location.  It  is  conjectured  that  forcing  at  the  higher 
frequency  kept  the  flow  more  periodic,  thus  reducing  the  vortex  pairing  tendency,  which  created 
the  largest  structures  and  therefore  the  largest  optical  distortions.  Interestingly,  open  loop  forcing 
at  these  higher  frequencies  did  not  show  this  level  of  performance,  which  was  attributed  to  phase 
and  frequency  differences  between  flow  states  and  forcing  input.  Only  with  feedback  control  was 
it  possibly  to  react  to  these  differences  in  the  adaptive  manner  necessary  to  reduce  the  density 
fluctuations  in  the  shear  layer. 
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Figure  52:  Mode  amplitudes  and  control  output  of  CFD  closed  loop  simulation.  The  control 
output  is  shown  for  the  complete  simulation,  periodic  forcing  for  t  <  0.025s,  feedback  control  for 
0.025s  <  t  <  0.06s,  and  unforced  for  t  >  0.06s. 
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Figure  53:  OPD  calculation  of  closed  loop  CFD  simulation  with  adjusted  controller.  Periodic 
forcing  for  0s  <  t  <  0.02551,  closed  loop  simulation  for  0.025  <  t  <  0.06s,  and  unforced  for  t  > 
0.065.  Reduction  of  OPD  on  the  order  of  30%  of  the  OPD  is  seen. 
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5.4  Axisymmetric  Bluff  Body 

5.4.1  Overview 

Classical  control  theory  is  very  limited  when  dealing  with  high-dimensional,  extremely  non-linear 
systems  such  as  flow  fields.  Flow  fields  are  governed  by  the  Navier  Stokes  equations,  Eq  74,  a  set 
of  second  order,  non-linear  partial  differential  equations.  New  techniques  need  to  be  established  to 
make  use  of  current  control  theories,  while  also  allowing  for  a  reasonable  design  process  for  linear, 
non-linear,  or  adaptive  control  for  complex  flow  fields. 

du 

dt 

The  synopsis  of  active  feedback  flow  control  is  to  use  a  fluidic  actuator  on  an  aerodynamics 
body  which  is  able  to  perturb  the  flow  away  from  the  original  state  and  typically  cause  some  type 
of  desired  response,  for  example  increased  lift  coefficient,  regulation  of  undersized  loads,  drag 
reduction,  optical  effects,  vortex  positioning,  etc.  Sensors  on  the  body  measure  the  flow  state 
which  is  then  translated  into  an  actuation  input  through  some  control  algorithm.  This  is  shown 
by  an  example  to  a  forebody  at  high  angle  of  attack  in  Figure  54.  The  challenges  associated  with 
active  feedback  flow  control  are  actuator  placement,  sensor  placement  and  model/control  algorithm 
design. 


+  u  •  Vu  =  —V p  +  V  •  T  +  f, 


(74) 


Filter 
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Mode 
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Figure  54:  Flow  control  approach  used  for  design  and  implementation  of  reduced  order  model 
based  control 

Figure  55  shows  the  approach  adopted  by  the  USAFA  flow  control  research  group.  This  frame¬ 
work  is  a  systematic  road  map  to  developing  a  reduced  order  model,  control  algorithm,  and  optimal 
sensor  placement  for  non-linear  fluid  dynamic  systems.  The  ultimate  goal  of  flow  control  research 
is  to  develop  a  robust  model  and  control  algorithm  for  a  specific  flow  field  to  provide  that  as  a 
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deliverable  to  the  customer.  The  secondary  goal  of  this  flow  control  research  is  to  gain  physical 
insight  into  the  fluid  dynamics  through  closing  the  loop. 


Figure  55:  Flow  control  approach  used  for  design  and  implementation  of  reduced  order  model 
based  control 

In  the  past  it  has  been  seen  that  the  more  insight  and  understanding  of  the  fluidic  mechanisms 
at  play  increase  the  fidelity  of  the  reduced  order  model  and  performance  of  the  control  algorithm. 
For  instance  while  closing  the  loop  on  the  backward  facing  step,  the  closed  loop  dynamics  showed 
that  exciting  the  Kelvin  Helmholtz  structures  in  a  certain  location  actually  caused  them  to  dissipate 
just  afterwards.  Thus,  the  optical  abberations  were  able  to  reduced  over  an  aperture. 

The  approach  begins  with  developing  a  fluidic  actuator.  Types  of  actuators  consist  of  SDBD 
plasma,  blowing  and  suction  ports,  synthetic  jets,  flow  effectors,  speakers,  etc.  Placement  of  actu¬ 
ators  and  number  of  actuators  is  typically  chosen  by  rules  of  thumb  from  fluid  dynamics.  This  is 
a  suboptimal  approach.  For  future  success  of  feedback  flow  control,  an  autonomous  method  needs 
to  be  determined  for  actuator  placement  and  design.  This  program  will  evaluate  different  technical 
solutions  to  this  issue  on  the  forebody  flow  control  problem. 

Once  the  fluidic  actuator  is  in  place  either  in  CFD  or  experiments,  open  loop  dynamics  are 
acquired  through  various  forms  of  input.  Typically,  step,  impulse,  ramp,  and  periodic  inputs  are 
used  to  quantify  the  system  dynamics.  As  shown  in  previous  applications  of  feedback  flow  control, 
these  open  loop  dynamics  are  a  crucial  step  in  understanding  the  flow  field.  The  actuator  to  fluidic 
response,  i.e.  the  controllability  of  the  fluidic  system,  is  a  critical  relationship  which  is  essential. 
For  unsteady  flow  fields,  the  state  trajectories  from  un-actuated  to  actuated  states  and  vice  versa 
is  non-unique  and  highly  dependent  upon  the  initial  state  of  the  flow.  This  necessitates  an  optimal 


82 


Fagley 


KC  Engineering,  Inc. 


April  8,  2013 


FA7000- 1 0-2-0003 .Final 


selection  of  open-loop  forcing  which  pseudo-randomly  captures  each  state  trajectory.  This  infor¬ 
mation  will  then  provide  for  the  compilation  of  a  flow  state  database,  which  is  essential  for  the 
development  of  a  reduced  order  model  which  accurately  represents  the  unforced,  open-loop,  and 
transient  states  near  the  desired  controlled  state.  While  these  questions  are  best  answered  using 
CFD  simulations  because  of  the  detailed  data  available,  it  is  impossible  to  interrogate  the  whole 
parameter  space  this  way.  Experimental  investigations  are  also  necessary  to  provide  the  crucial 
survey  of  the  parameter  space  to  highlight  regions  of  particular  interest. 

Once  the  open  loop  database  is  formulated,  numerical  decomposition  techniques  are  used  to  ex¬ 
tract  pertinent  dynamics.  Specifically  these  decompositions  decouple  spatial  and  temporal  modes 
in  an  optimal  fashion.  For  the  development  of  Reduced  Order  Models  (ROMs)  of  the  flow  field,  a 
software  suite  developed  in  the  US  Air  Force  Department  of  Aeronautics  by  the  researcher  will  be 
readily  available  and  applied  to  new  problems.  The  data  analysis  part  of  the  software  suite  consists 
of  many  tools  such  as:  proper  orthogonal  decomposition  (POD),  double  POD,  balanced  POD,  dy¬ 
namic  modal  decomposition  (also  referred  to  as  Koopman  analysis)  and  wavelet  decompositions. 

Each  decomposition  has  unique  advantages  and  disadvantages  while  the  overarching  goal  is 
the  same  -  to  extract  the  dynamical  behavior  of  large  scale,  coherent  structures  in  the  flow  while 
decoupling  spatial  and  temporal  information.  Extracting  the  desired  or  dominant  dynamics  of  a 
fluid  field  is  a  highly  debated  topic.  For  instance  POD  defines  the  dynamics  through  largest  ener¬ 
getic  modes;  BPOD  defines  the  dynamical  modes  as  a  set  which  optimizes  the  observability  and 
controllability  grammians;  DPOD  emphasizes  an  energetic  mode  set  coupled  with  shift  or  pertur¬ 
bation  modes  to  the  dominant  modes  due  to  actuated  and  un-actuated  transients,  and  finally  the 
Koopman  analysis  extracts  spatial  growth  and  decay  rates  (globally  unstable  and  stable  modes)  as 
well  as  spatial  frequencies  (marginally  stable  modes).  Each  of  these  tools,  while  vastly  different 
mathematical  procedures  produce  the  same  result,  decoupled  spatial  and  temporal  information. 
These  decomposition  techniques  are  commonly  understood  for  unsteady  flow  fields,  but  the  exten¬ 
sion  and  application  to  deformable  bodies  in  computational  simulation  or  experimental  testing  has 
never  been  attempted.  It  is  the  intent  of  this  research  program  to  find  a  suitable  strategy  or  combi¬ 
nation  of  strategies  to  extract  the  dynamics  of  the  fluid,  structure,  and  fluidic  actuation  interaction. 

Once  an  understanding  of  the  underlying  physics  of  fluid  structure  interactions  is  produced, 
a  proper  dynamic  reduced  order  model  can  be  developed.  With  this  accurate  system  model  the 
optimal  type  and  distribution  of  actuators  and  surface  sensors  can  be  determined  and  implemented 
in  both  the  simulations  and  experiments.  During  previous  research  projects,  the  flow  control  group 
at  USAFA  has  found  that  the  low-dimensional  modeling  approach  is  the  most  beneficial  when  it 
comes  to  realizing  a  structured  model-based  closed-loop  strategy  for  flow  control.  Assuming  a 
suitable  mode  set  is  determined  from  the  previous  section  which  represents  the  unforced,  open- 
loop,  and  transient  states,  the  associated  temporal  dynamics  need  to  be  modeled. 

A  widely  accepted  approach  to  model  experimental  and  computational  periodic  flows  is  the 
Galerkin  Method.  A  Galerkin  projection  is  a  method  for  obtaining  approximations  to  a  high  dimen¬ 
sional  dynamical  system  by  projecting  the  underlying  dynamics  onto  a  reduced  order  subspace.  In 
the  application  to  fluid  dynamics,  the  Navier  Stokes  equations  are  projected  onto  a  subspace  which 
is  spanned  by  an  orthogonal  basis  which  represents  a  majority  of  the  system  dynamics  typically 
obtained  by  POD.  The  resulting  equations  are  a  set  of  non-linear  ODE’s  in  the  form  of, 

Xf  =  fo  +  Lxf(t)  +  Q(xf(t)(^}Xf(t)),  (75) 

where  the  linear  term  is  representative  of  the  viscous  term  in  the  Navier  Stokes  and  the  quadratic 


83 


Fagley 


KC  Engineering,  Inc. 


April  8,  2013 


FA7000- 1 0-2-0003 .Final 


term  is  representative  of  the  convective  and  pressure  terms  in  the  Navier  Stokes.  In  some  instances 
it  is  beneficial  to  also  add  a  higher  order  cubic  term  to  account  for  mean  flow  perturbations  and 
experimental  noise.  An  extension  of  the  Galerkin  Method  is  the  Discontinuous  Galerkin  Method 
which  may  be  essential  for  producing  a  stable  discretization  of  the  convective  operator  over  un¬ 
structured,  deformable  meshes. 

Previously,  this  research  investigated  wavelet  basis  networks  (WNARX)  to  demonstrate  the  full 
capability  of  identifying  complex  flow  response  for  a  range  of  open  loop  parameters.  The  WNARX 
represents  a  dynamic  model  which  can  simulate  off  design  flow  cases,  serve  as  reference  signal, 
and  ultimately  predict  closed  loop  behavior  for  control  design.  The  WNARX  model  uses  the  same 
network  architecture  as  a  neural  network;  the  only  difference  is  radial  basis  functions  are  used  as 
each  neuron’s  transfer  function.  This  is  shown  by,  'F(f)„i,s  =  *P  {—f),  where  u  is  the  translation  of 
the  wavelet,  s  is  the  dilation,  and  *P  is  referred  to  as  the  mother  wavelet,  which  is  a  radial  basis 
function  in  this  case.  WNARX  models  are  much  better  suited  for  identifying  the  frequency  rich 
dynamics  of  complex,  turbulent  flow  fields.  The  overall  WNARX  model  is  given  by, 

N 

fit)  =  J2w;-vP(5/(r-M;-))+c7’t+/o  (76) 

1=1 

where  Wj  are  the  weights,  N  is  the  number  of  wavelet  functions,  cT  represent  the  linear  connections, 
and  /o  is  the  bias.  This  proposal  will  use  this  new  system  identification  method  to  formulate  an 
extremely  low  dimensional  model  based  from  CFD  simulations  and  POD/DPOD  decompositions 
to  accurately  predict  closed  loop  dynamics  of  a  given  flow  field.  This  model  is  then  used  to  perform 
feedback  simulations  to  condition  control  algorithms  which  can  then  be  applied  directly  to  CFD 
simulations  and  experiments. 

Once  this  model  (either  Galerkin  model,  or  WNARX)  is  validated  over  unforced  and  open- 
loop  parameter  space,  it  will  serve  as  the  basis  for  the  control  strategy  development.  Typical 
control  approaches  to  non-linear  systems  can  be  used.  For  instance  direct  adaptive  control,  sliding 
mode  control,  or  model  predictive  control  may  be  used  to  achieve  stable,  robust  reference  tracking 
ability.  Previous  research  supports  that  adaptive  control  adequately  handles  model  uncertainties 
between  model  and  CFD  closed  loop  simulations.  This  research  project  aims  at  evaluating  the  per¬ 
formance  and  stability  criterion  of  these  three  control  approaches.  As  for  differences  between  CFD 
simulations  and  experimental  tests,  the  actuation  dynamics  and  sensor  dynamics  will  be  unique  to 
each  environment,  but  the  underlying  control  theory,  flow/structure  state  definition,  and  estima¬ 
tion  algorithms  will  remain  constant.  This  will  allow  for  model  and  estimation  development  to 
be  based  upon  both  numerical  and  experimental  data.  In  a  sense  the  control  will  be  modularized 
by  the  actuation/sensing  dynamics  which  will  provide  matching  CFD  and  experimental  results  in 
closed  loop  simulation/experimentation. 
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5.4.2  Open  Loop  Dynamics  Experiment 

As  presented  in  Farnsworth  et  al.  [2012],  at  the  current  operating  conditions,  the  system  is  shown 
to  be  proportional  at  an  incidence  of  50°.  Figure  56  shows  a  schematic  of  the  responsiveness 
of  the  asymmetric  state  (measured  as  the  resulting  side  force)  to  plasma  actuation.  The  positive 
x-axis  represents  the  port  forcing  strength  and  the  negative  x-axis  represents  the  starboard  forc¬ 
ing  strength;  the  zero  location  is  the  unforced  state.  The  unforced  state  of  the  asymmetric  vortex 
configuration  varies  based  on  geometry  disturbances,  flow  conditions,  misalignments,  flow  imper¬ 
fections  etc.  Around  this  initial  state  is  a  dead  zone  in  the  actuator  dynamics;  that  is,  the  actuation 
voltage  must  exceed  a  certain  limit  before  plasma  formation  takes  place.  Above  and  below  this  re¬ 
gion  a  linear  response  in  asymmetric  vortex  state  was  found.  At  large  enough  forcing  magnitudes 
the  vortex  system  does  saturate  in  the  fully  left  or  right  asymmetric  vortex  state. 


Figure  56:  Representation  of  the  forcing  characteristics  of  the  asymmetric  vortex  state  due  to 
plasma  actuation. 

Figure  56  is  experimentally  verified  by  side  force  and  sectional  pressure  measurements  as  de¬ 
picted  in  Figure  57.  The  side  force,  Cy  and  sectional  pressure  coefficient,  A Cp,  at  x/D  =  3  vary 
analogously  with  varying  port/starboard  plasma  voltage  which  supports  that  time  resolved  pres¬ 
sure  measurements  do  accurately  correlate  with  integrated  force  measurements.  As  Figure  57  also 
shows,  the  system  responds  nearly  proportionally,  although  non-linear  effects  are  apparent.  For 
instance,  the  dead  zone  in  the  actuation  voltage  range  from  —5 kV  <V<  5 kV  does  exist;  this  is, 
primarily  due  to  the  fact  that  plasma  has  not  formed  at  these  smaller  voltage  potentials.  Also, 
a  hysteresis  is  definitely  observed,  that  is  the  path  along  a  positive  voltage  gradient,  W  >  °’  is 
different  from  the  path  along  a  negative  voltage  gradient,  ^  <  0,  in  both  time  accurate  (red/blue 
lines)  and  integrated  measurements  (black  lines).  Finally,  a  larger  gradient,  (-^L  is  seen  near  the 
symmetric  vortex  location  ( Cy  =  A  Cp  =  0)  ,  supporting  the  fact  that  a  small  amount  of  bistability 
does  exist  at  this  flow  condition. 
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Figure  57:  (a)  Side  force  coefficient  and  (b)  differential  pressure  coefficient  at  x/D  =  3  versus  the 
plasma  actuator  voltage  during  the  ramp  modulation. 

The  system  response  in  Figure  57  is  well  modeled  by  a  linear  system  with  slight  non-linearities 
present. 

5.4.2.1  Step  Response  The  response  of  the  asymmetric  vortex  state  due  to  a  plasma  actua¬ 
tion  step  input  is  used  to  develop  a  linear  time  invariant  model.  The  step  response  test  campaign 
consisted  of  a  modulated  square  wave  at  a  frequency  of  1  Hz  for  a  total  of  20  periods.  The  data 
was  then  phase  averaged  over  a  test  duration  of  20  seconds  to  reduce  measurement  noise.  The 
amplitude  of  the  step  was  at  maximum  operational  voltage  of  12  kV  before  the  amplifier  began 
displaying  non-linear  effects.  Figure  58  shows  the  normalized  response  to  the  step  input  at  initial 
transient  times  and  Figure  59  for  the  ending  transient  times. 

The  overall  time  delay  consists  of  the  convective  time  delay  for  the  disturbance  to  reach  the 
sensor  location,  lag  time  for  the  fluid  to  respond,  and  transition  time  to  achieve  90%  of  the  steady- 
state  value.  To  decouple  each  of  these  sources  of  delay,  the  convective  time  delay  is  defined  as 
the  time  from  which  the  step  begins  to  the  time  at  which  a  10%  change  in  the  unforced  steady- 
state  value  is  observed.  The  rise  time  is  defined  as  the  time  from  a  10%  change  in  the  unforced 
steady  state  value  to  the  time  at  which  90%  of  the  forced  steady-state  value  was  achieved.  The 
time  responses  are  then  normalized  by  the  flow  through  time,  t,  which  is  defined  as, 


L, 


' cone 


(77) 


and  was  measured  to  be  approximately  13  ms  at  the  current  operating  conditions.  These  times 
are  summarized  for  the  rising  and  falling  transients  in  Table  7.  The  lag  time  or  presence  of  non¬ 
minimum  phase  are  difficult  issues  to  decouple  in  the  dynamics  so  further  analysis  techniques  are 
necessary. 
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Figure  58:  Initial  transient  of  the  vortex  state  due  to  step  input.  Linear  combination  of  pressure 
measurements  at  x/D  =  2  and  x/D  =  3  to  estimate  the  side  force  shown  in  green.  Blue  shows  the 
step  change  of  the  actuation  input  in  kV.  Also,  the  convective  delay  time,  Trci,  is  shown  in  cyan  and 
the  transition/rise  time,  Trt,  is  shown  in  red. 


Figure  59:  Ending  transient  of  the  vortex  state  due  to  step  input.  Linear  combination  of  pressure 
measurements  at  x/D  =  2  and  x/D  =  3  to  estimate  the  side  force  shown  in  green.  Blue  shows  the 
step  change  of  the  actuation  input  in  kV.  Also,  the  convective  delay  time,  Tfj,  is  shown  in  cyan  and 
the  transition/rise  time,  Tft,  is  shown  in  red. 


Table  7:  Rise  and  fall  time  summary 


Initial  Transient  (Port  -  Starboard)  Time  [?*] 

Ending  Transient  (Starboard-Port)  Time  [ t *] 

Delay  Time 

Rise  Time 

Total 

%  Overshoot 

Delay  Time 

Rise  Time 

Total 

%  Overshoot 

Cy 

0.6375 

0.69 

1.3275 

18.5  % 

0.6070 

0.57 

1.177 

24.4% 
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As  shown  in  Table  7,  the  convective  times  for  a  port  to  starboard  or  vice  versa  are  very  consis¬ 
tent,  as  expected.  The  rise  or  transition  time  from  starboard  to  port  is  a  bit  faster  than  the  transition 
time  from  port  to  starboard.  This  is,  arguably,  because  the  initial  asymmetric  state  prefers  the  port 
side  due  to  geometry  imperfections,  flow  misalignments,  etc.;  thus  the  flow  prefers  transition  back 
to  the  port  state  and  induces  a  restoring  force,  reducing  transient  time  and  also  causing  a  larger 
overshoot  of  the  steady  state  value.  Nonetheless,  the  dynamics  of  the  asymmetric  vortex  problem 
as  shown  by  the  step  response  are  very  well  represented  by  a  linear  time  invariant  system. 

5.4.2.2  Sinusoidal  Response  To  determine  the  frequency  response  of  the  system  dynamics, 
sinusoidal  forcing  is  used  on  the  port  actuator  and  the  system  response  is  observed  by  the  linear 
combination  of  all  of  pressure  measurements  as  given  by  Cy.  The  actuation  voltage  is  modulated 
by  an  offset  sinusoid,  by  the  equation, 

A(t)  =  Vmax(siii(2,7r  (Qt^  +0.5). 

The  test  durations  consisted  of  30  seconds  with  a  sampling  frequency  of  10  kHz.  An  example  of 
forcing  at  a  frequency  of  20  Hz  is  shown  in  Figure  60.  The  input-output  signals  are  shown  as  well 
as  the  frequency  spectrum.  Figure  60b  shows  a  large  peak  at  the  forcing  frequency  showing  the 
fluidic  receptivity  to  the  forcing.  A  frequency  sweep  was  conducted  over  a  wide  frequency  range 
to  determine  the  cutoff  frequency  as  well  as  the  magnitude  and  phase  of  the  system. 

The  natural  rise  time  of  the  fluidic  response  is  approximately  1.1 1*  1 .6+,  depending  upon 
port  to  starboard  actuation  or  vice  versa,  due  to  a  unit  step  input  as  shown  in  section  Table  7. 
The  natural  frequency  is  approximately  50  Hz.  This  suggests  that  a  pole  exists  near  this  location. 
Because  of  this  observation,  the  modulation  frequency  was  chosen  at  discrete  locations  over  the 
range  of  0.1  Hz  <  ft)  <  200Hz,  to  determine  magnitude,  phase  and  cutoff  frequency. 


Forcing,  F  =  0.27  1/t 


Figure  60:  Pressure  location  x/D  =  2  (a)  Time  domain  forcing  and  response  data  for  Frequency  = 
20  Hz  (b)  Time  domain  forcing  and  response  data  for  Frequency  =  20  Hz. 

For  all  of  the  forcing  frequencies  the  data  is  summarized  in  Figure  88  where  the  experimental 
data  is  plotted  in  red.  From  this  frequency  response  information  the  cutoff  frequency  can  be  esti¬ 
mated  by  a  -3  dB  attenuation  point.  This  is  computed  to  be  approximately  50  Hz  and  corresponds 
to  a  80  degree  phase  lag. 
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Figure  61:  Experimental  frequency  and  phase  response  information  for  sinusoidal  forcing  cam¬ 
paign 

5.4.2.3  Impulse  Response  The  impulse  response  of  the  asymmetric  vortex  state  was  also  mea¬ 
sured.  For  these  open-loop  tests  the  duty  cycle  was  varied  for  a  square  modulation  wave  at  a 
frequency  of  10 Hz  over  a  range  of  1%  to  20%.  The  experimental  measurements  are  shown  below 
in  Figure  89  for  the  different  duty  cycles.  The  initial  flow  state  was  shifted  to  Cy  =  0,  i.e.  the 
symmetric  state  for  modeling  purposes.  All  of  the  impulses  were  initiated  at  t  =  0,  so  that  the 
flow  response  is  aligned  for  each  duty  cycle.  These  measurements  were  phase  averaged  over  a 
100  total  cycles.  The  results  in  Figure  89  are  well  depicted  by  a  linear  system.  As  the  duty  cycle 
increases  beyond  10%  an  amount  of  undershoot  is  seen  by  the  vortex  dynamics.  This  data  set  is 
used  completely  as  validation  for  the  model  formulation  and  model  selection  technique. 
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Figure  62:  Experimental  measurements  of  the  impulse  response  with  varying  duty  cycle  of  the 
pulse  width  at  a  maximum  voltage  of  1 .4  kV.  This  data  serves  purely  as  validation  data  for  model 
development  in  subsequent  sections. 

5.4.3  Open  Loop  Dynamics  Simulations 

As  presented  in  Farnsworth  et  al.  [2012],  at  the  current  operating  conditions,  the  system  is  shown 
to  be  proportional  at  an  incidence  of  50°.  Figure  63  shows  a  schematic  of  the  responsiveness 
of  the  asymmetric  state  (measured  as  the  resulting  side  force)  to  plasma  actuation.  The  positive 
x-axis  represents  the  port  forcing  strength  and  the  negative  x-axis  represents  the  starboard  forc¬ 
ing  strength;  the  zero  location  is  the  unforced  state.  The  unforced  state  of  the  asymmetric  vortex 
configuration  varies  based  on  geometry  disturbances,  flow  conditions,  misalignments,  flow  imper¬ 
fections  etc.  Around  this  initial  state  is  a  dead  zone  in  the  actuator  dynamics;  that  is,  the  actuation 
voltage  must  exceed  a  certain  limit  before  plasma  formation  takes  place.  Above  and  below  this  re¬ 
gion  a  linear  response  in  asymmetric  vortex  state  was  found.  At  large  enough  forcing  magnitudes 
the  vortex  system  does  saturate  in  the  fully  left  or  right  asymmetric  vortex  state. 
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Figure  63:  Representation  of  the  forcing  characteristics  of  the  asymmetric  vortex  state  due  to 
plasma  actuation. 

Figure  63  is  verified  in  CFD  simulations  by  side  force  as  depicted  in  Figure  64.  The  side  force, 
Cy  varies  proportionally  with  varying  port/starboard  plasma  voltage.  As  Figure  64  also  shows,  the 
system  responds  nearly  proportionally,  although  non-linear  effects  are  apparent.  For  instance,  the 
dead  zone  in  the  actuation  strength  from  — 0.025g/s  <m<  0.025g/s  does  exist.  Saturation  regions 
are  also  seen  for  m  >  0.75 g/s. 

The  system  response  in  Figure  64  is  well  modeled  by  a  linear  system  with  slight  non-linearities 
present. 


5.4.3.1  Step  Response  The  response  of  the  asymmetric  vortex  state  due  to  a  plasma  actua¬ 
tion  step  input  is  used  to  develop  a  linear  time  invariant  model.  The  step  response  test  campaign 
consisted  of  a  modulated  square  wave  at  a  frequency  of  1  Hz  for  a  total  of  20  periods.  The  data 
was  then  phase  averaged  over  a  test  duration  of  20  seconds  to  reduce  measurement  noise.  The 
amplitude  of  the  step  was  at  maximum  operational  voltage  of  12  kV  before  the  amplifier  began 
displaying  non-linear  effects.  Figure  65  shows  the  normalized  response  to  the  step  input  at  initial 
transient  times  and  Figure  66  for  the  ending  transient  times. 

The  overall  time  delay  consists  of  the  convective  time  delay  for  the  disturbance  to  reach  the 
sensor  location,  lag  time  for  the  fluid  to  respond,  and  transition  time  to  achieve  90%  of  the  steady- 
state  value.  To  decouple  each  of  these  sources  of  delay,  the  convective  time  delay  is  defined  as 
the  time  from  which  the  step  begins  to  the  time  at  which  a  10%  change  in  the  unforced  steady- 
state  value  is  observed.  The  rise  time  is  defined  as  the  time  from  a  10%  change  in  the  unforced 
steady  state  value  to  the  time  at  which  90%  of  the  forced  steady-state  value  was  achieved.  The 
time  responses  are  then  normalized  by  the  flow  through  time,  T,  which  is  defined  as, 


T  = 


Leone 

-f/oT’ 


(78) 


and  was  measured  to  be  approximately  13  ms  at  the  current  operating  conditions.  These  times 
are  summarized  for  the  rising  and  falling  transients  in  Table  8.  The  lag  time  or  presence  of  non- 
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minimum  phase  are  difficult  issues  to  decouple  in  the  dynamics  so  further  analysis  techniques  are 
necessary. 


Table  8:  Rise  time  summary 


Initial  Transient  (Port  - 

Starboard)  Time  [t] 

Delay  Time 

Rise  Time 

Total 

%  Overshoot 

A 

Cy 

0.62 

0.39 

1.01 

15.5  % 

When  comparing  CFD  transient  times  to  experimental  transient  times,  very  good  agreement 
exists,  see  Table  8  and  Table  7.  The  simulations  are  slightly  faster  than  the  experiments  (  30%). 
This  is  mainly  due  to  transient  dynamics  of  the  experimental  circuitry  which  is  non  existent  in  the 
CFD  simulation.  Nonetheless,  the  dynamics  of  the  asymmetric  vortex  problem  as  shown  by  the 
step  response  are  very  well  represented  by  a  linear  time  invariant  system. 

5.4.4  Impulse  Response 

The  impulse  response  of  the  asymmetric  vortex  state  was  also  simulated.  For  these  open-loop  tests 
the  duty  cycle  was  varied  for  a  square  modulation  wave  at  a  frequency  of  10 Hz  over  a  range  of 
1%  to  10%.  The  simulation  results  are  shown  below  in  Figure  67  for  the  different  duty  cycles. 
The  initial  flow  state  was  shifted  to  Cy  =  0,  i.e.  the  symmetric  state  for  modeling  purposes.  All 
of  the  impulses  were  initiated  at  t  =  0,  so  that  the  flow  response  is  aligned  for  each  duty  cycle. 
The  results  in  Figure  67  are  well  depicted  by  a  linear  system.  This  data  set  is  used  completely  as 
validation  for  the  model  formulation  and  model  selection  technique. 


Figure  67:  Impulse  response  of  the  CFD  simulation  varying  duty  cycle  of  the  pulse  width  at  a  unity 
magnitude  of  C ^ .  This  data  serves  purely  as  validation  data  for  model  development  in  subsequent 
sections. 
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Figure  64:  Average  side  force  coefficient  for  different  forcing  parameters.  Port  forcing  corresponds 
to  a  negative  x-values  while  starboard  forcing  corresponds  to  a  positive  x-values. 
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Figure  65:  Side  force  in  green  shows  the  response  due  to  a  step  input  in  C),  in  blue. 


Figure  66:  Initial  transient  of  the  vortex  state  due  to  step  input.  Side  force  in  green  shows  the 
response  due  to  a  step  input  in  C ^  in  blue.  Also,  the  convective  delay  time,  Tfd,  is  shown  in  cyan 
and  the  transition/rise  time,  Tft,  is  shown  in  red. 


94 


Fagley 


KC  Engineering,  Inc. 


April  8,  2013 


FA7000- 1 0-2-0003 .Final 


5.4.5  Numerical  decomposition  and  flow  state  definition 

To  further  investigate  the  flow  field  around  a  von  Karman  tangent  ogive  and  develop  a  Reduced- 
Order-Model  (ROM)  for  feedback  flow-control,  unsteady  numerical  investigations  were  under¬ 
taken.  The  simulations  are  performed  using  Cobalt,  an  unstructured  finite- volume  code  developed 
for  the  solution  of  the  compressible  Navier-Stokes  equations.  The  basic  algorithm  is  described  in 
Strang  et  al., Strang  et  al.  [1999]  although  substantial  improvements  have  been  made  since  then. 
The  numerical  method  is  a  cell-centered  finite  volume  approach  applicable  to  arbitrary  cell  topolo¬ 
gies  (e.g,  hexahedra,  prisms,  tetrahedra).  The  spatial  operator  uses  a  Riemann  solver,  least  squares 
gradient  calculation  with  QR  factorization  to  provide  second  order  accuracy  in  space.  A  point 
implicit  method  using  analytic  first-order  inviscid  and  viscous  Jacobians  is  used  for  advancement 
of  the  discretized  system.  For  time-accurate  computations,  a  Newton  sub-iteration  scheme  is  em¬ 
ployed,  resulting  in  a  method  that  is  formally  second  order  accurate  in  time.  For  parallel  per¬ 
formance,  Cobalt  utilizes  the  domain  decomposition  library  ParMETIS  to  provide  optimal  load 
balancing  with  a  minimal  interface  between  zones. Karypis  et  al.  [1997] 

5.4.5.1  Grid  and  Model  Geometry  The  geometry  considered  in  this  investigation  is  a  generic 
tangent  ogive  with  fineness  ratio  fr  =  3.5  and  a  model  base  diameter  of  D  =  0.1  m.  At  the  base 
of  the  ogive,  a  0.05  m  long  cylindrical  section  has  been  added  such  that  the  overall  length  of 
the  model  is  0.40  m.  This  model  geometry  was  chosen  to  match  an  accompanying  wind  tunnel 
experiment. Fagley  et  al.  [2012a]  For  reference,  the  center  of  the  coordinate  system  is  at  the  nose  of 
the  model.  The  positive  x-direction  extends  along  the  body,  the  positive  y-direction  points  in  the 
starboard  direction,  and  the  positive  z-direction  is  up,  normal  to  the  body.  For  the  simulations,  all 
the  reference  conditions  are  set  to  standard  sea  level,  with  an  inflow  Mach  number  of  M=0.1;  this 
results  in  a  Reynolds  number  based  on  the  base  diameter  of  Re  =  220, 000. 

The  grid  was  generated  using  Simcenter/SolidMesh  and  had  approximately  16M  cells  (Fig.  68). 
To  avoid  asymmetries  in  the  flow  field  as  a  result  of  an  asymmetric  grid,  the  grid  was  generated 
around  half  the  model  and  then  mirrored;  therefore,  the  grid  is  symmetric  on  the  port  and  starboard 
sides  of  the  model.  Two  patches  to  simulate  the  plasma  actuators  used  in  the  accompanying  exper¬ 
iment  were  added  to  the  model  at  90°  from  the  model’s  meridian.  The  start  of  the  boundary  patch 
was  placed  0.4  cm  from  the  tip  of  the  model  and  was  1  cm  long.  All  model  boundaries  were  always 
set  to  solid-wall,  no-slip  conditions  except  during  the  open-loop  forcing  investigations  where  one 
of  the  actuator  boundary  patch  (either  port  or  starboard)  was  switched  to  a  moving  wall  bound¬ 
ary  patch.  In  the  open-loop  simulations,  the  moving  wall  boundary  patch  prescribed  a  tangential 
velocity  at  the  wall  in  the  direction  of  the  freestream.  A  spherical  farfield  boundary  was  placed 
40  diameters  away  from  the  model  to  minimize  the  influence  of  pressure  reflections.  For  all  the 
calculations,  the  farfield  used  a  modified  Riemann  condition  and  the  time  step  was  specified  at 
At  =  0.0001  s. 
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Figure  68:  Ogive  geometry  and  grid  with  simulated  port  actuator  shown  in  purple:  a)  front  view 
with  grid  shown  at  the  center  plane  of  the  forcing  patch,  b)  side  view  of  the  grid  along  the  centerline 

the  time-resolved  simulation  data  around  the  ogive  body,  an  array  of  Cobalt  ’taps’ 
was  used.  Figure  69  shows  the  tap  grid  used  around  the  ogive  model.  The  tap  grid  extended  along 
the  whole  body,  from  x/D  =  0  to  x/D  =  4  and  110°  from  the  leeward  meridian.  The  taps  were 
non-uniformly  spaced,  such  that  a  high  spatial  resolution  is  obtained  near  the  body,  especially  near 
the  tip  of  the  model.  In  total,  37,625  taps  were  used  to  extract  all  the  flow  quantities  during  each 
time  iteration. 


(a)  (b) 

Figure  69:  Ogive  geometry  with  tap  grid:  a)  isometric  view  and  b)  front  view. 


5.4.5.2  Proper  Orthogonal  Decomposition  For  data  analysis,  Proper  Orthogonal  Decompo¬ 
sition  (POD)  has  been  shown  to  be  a  very  effective  tool  to  characterize  flow  fields. Sirovich  [1987], 
Berkooz  et  al.  [1993]  In  POD,  highly  complex  flow  fields  are  decomposed  into  spatial  modes  with 
a  corresponding  time  varying  amplitude  (or  coefficients): 

K 

(p(x,y,t )  =  £  ak(t)<j)k(x,y)  (79) 

k=l 
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where  (p  is  a  flow  quantity,  <j>k{x,y)  are  the  spatial  modes,  and  ak(t)  are  the  time  coefficients. 
While  this  decomposition  is  well  suited  to  time  periodic  flow  fields,  it  faces  problems  for  transient 
flowsSiegel  et  al.  [2005b]  and  potentially  aperiodic  flows.  In  these  cases,  POD  is  based  on  all 
the  snapshots  used  to  generate  the  modes,  minimizing  the  overall  error  with  the  fewest  number  of 
modes  possible.  Therefore,  transients  or  aperiodicity  in  the  flowfield  can  be  missed  due  to  the  small 
contribution  made  to  the  overall  estimates  of  the  spatial  and  temporal  modes.  Different  additions 
to  the  basic  POD  procedure  have  been  proposed,  most  notably  the  addition  of  a  shift  mode  as 
introduced  independently  by  Noack  Noack  et  al.  [2003]  as  well  as  Siegel  et  al.Siegel  et  al.  [2003b] 

5.4.5.3  Results  Before  investigating  the  effects  of  open-loop  forcing  on  the  side  force  gener¬ 
ated  from  the  asymmetric  vortex  state,  unforced  simulations  were  carried  out  where  both  the  port 
and  starboard  forcing  patches  were  set  to  solid-wall,  no-slip  boundary  conditions.  In  this  config¬ 
uration,  the  model  is  geometrically  perfect  such  that  no  disturbance  exists  to  cause  the  leeward 
vortices  to  lock  into  an  asymmetric  state.  Figure  70  shows  the  resultant  side  force  on  the  model  as 
a  function  of  time.  As  shown,  the  side  force  coefficient,  Cy,  fluctuates  between  the  port  (negative 
Cy  values)  and  starboard  (positive  Cy  values)  direction  reaching  magnitudes  slightly  greater  than 
0.5.  As  summarized  by  Bridges  et  al., Bridges  [2006]  if  the  root  cause  of  the  asymmetric  vortex 
configuration  is  the  result  of  a  convective  instability,  numerical  codes  should  not  produce  an  asym¬ 
metric  wake  on  geometrically  perfect  bodies.  If  they  do,  then  the  potential  exist  that  the  asymmetry 
is  due  to  numerical  issues  and  is  not  necessarily  real;  in  other  words,  the  right  solution  is  obtained 
for  the  wrong  reasons.  Figure  70  shows  that  the  geometrically  perfect  model  produces  an  average 
side  force  of  zero  using  the  current  numerical  setup.  To  verify  that  the  numerical  code  would  pro¬ 
duce  an  asymmetric  wake  when  a  geometrical  disturbance  was  present,  simulations  were  also  run 
with  a  small  bump  (0.5  mm  diameter  pin,  0.5  mm  tall)  placed  on  the  port  side  near  the  tip  of  the 
model.  Figure  70  also  shows  the  resultant  side  force  on  the  model  with  this  geometric  disturbance 
present.  As  shown,  the  geometric  disturbance  causes  the  side  force  on  the  model  to  lock  into  one 
side  (in  the  opposite  direction  of  the  disturbance)  and  the  magnitude  of  the  asymmetry  increases. 
Therefore,  at  this  Reynolds  number,  a  disturbance  on  the  port  side  causes  the  port  side  vortex  to 
lift  off  the  body,  causing  the  side  force  in  the  starboard  direction. 


Figure  70:  Comparison  of  the  side  force  coefficient  on  a  tangent  ogive  forebody  with  and  without 
a  geometric  disturbance  located  near  the  tip  of  the  model. 
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Following  the  unforced  simulations,  open-loop  simulations  which  modeled  the  plasma  actuator 
used  in  the  experiments  were  performed.  To  simulate  the  plasma  actuator,  a  moving  wall  boundary 
condition  was  implemented,  where  instead  of  setting  the  tangential  velocity  at  the  wall  to  zero  (i.e 
no  slip),  a  fixed  velocity  is  prescribed.  Therefore,  downstream  of  the  moving  wall  patch,  a  wall  jet 
is  formed  as  the  velocity  at  the  wall  diffuses  into  the  rest  of  the  boundary  layer.  The  prescribed  ve¬ 
locity  of  the  moving  wall  boundary  patch  was  changed  throughout  a  series  of  simulations,  ranging 
from  4  m/s  up  to  16  m/s  (10  to  45  percent  of  freestream).  While  these  velocities  are  much  greater 
then  those  generated  by  the  plasma  actuators  in  the  experiments,  the  resultant  wall  jets  downstream 
of  the  actuator  are  similar.Lee  et  al.  [2012] 

Figure  71  shows  the  resultant  side  force  on  the  model  for  a  moving  wall  velocity  of  12  m/s, 
where  the  actuation  is  switched  from  the  starboard  to  the  port  side  every  300  ms.  As  shown,  there 
is  a  delay  (~  100  ms)  between  the  time  actuation  is  turned  on  and  the  time  when  the  side  force  fully 
switches  sides  and  locks  into  its  new  value.  From  the  simulations,  it  takes  approximately  50  ms 
(five  convective  time  scalesLanser  and  Meyn  [1994])  before  the  side  force  begins  to  respond  to  the 
forcing  input.  The  other  50  ms  is  the  approximate  time  it  takes  for  the  vortices  on  the  model  to 
switch  states.  This  is  in  contrast  to  experimental  findings  based  on  dynamic  pressure  transducers 
on  the  model  which  show  that  the  model  responds  to  forcing  from  the  plasma  actuator  within  one 
to  two  convective  time  scales.Farnsworth  et  al.  [2012] 

Based  on  the  data  shown  in  Fig.  71,  the  side  on  which  forcing  is  applied  causes  the  vortex  on 
that  side  of  the  model  to  separate  from  the  body  and  the  side  force  is  directed  away  from  the  side  of 
actuation.  This  is  the  same  trend  observed  with  the  geometric  disturbance  introduced  on  the  model 
as  well  as  experimental  data  ?  for  a  similar  Reynolds  number  range.  This  indicates  that  at  this 
Reynolds  number,  instead  of  the  plasma  adding  momentum  to  the  flow  to  help  keep  the  boundary 
layer  attached,  it  is  instead  creating  a  disturbance  causing  the  flow  to  separate.  This  is  a  different 
trend  from  Matsuno  et  al.  Matsuno  et  al.  [2009]  who  tested  the  use  of  DBD  plasma  actuators  on 
a  ogive  model  at  a  lower  Reynolds  numbers  (~  50,000).  During  their  tests,  the  plasma  actuator 
caused  the  lifted  vortex  to  attach  to  the  body;  this  reattachment  of  the  lifted  vortex  was  attributed 
to  the  Coanda  effect.  At  the  lower  Reynolds  numbers  tested,  it  is  surmised  that  the  strength  of 
the  wall  jet  created  by  the  plasma  relative  to  the  freestream  is  much  larger  creating  a  Coanda  type 
effect,  while  at  the  higher  Reynolds  numbers,  the  jet  creates  a  disturbance  causing  the  flow  on  that 
side  to  separate. 
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Figure  71:  Open-loop  side  force  coefficient  on  a  tangent  ogive  body,  switching  between  starboard 
and  port  forcing  every  300  ms. 

To  illustrate  the  coherent  vortical  structures  in  the  flowfield,  Fig.  72  shows  an  isosurface  of 
Q  with  a  moving  wall  velocity  of  12  m/s  for  the  starboard  actuator  turned  on.Jeong  and  Hussain 
[1995]  As  shown,  at  this  Reynolds  number,  the  asymmetry  of  the  two  primary  vortices  is  small, 
even  though  the  side  force  coefficient  is  approximately  -1.0  at  this  instant.  To  help  distinguish 
between  the  port  and  starboard  vortices,  Fig.  72  is  colored  by  the  x-vorticity.  Looking  at  the  side 
view  of  Fig.  72,  smaller  coherent  structures  can  be  seen  feeding  into  the  primary  vortices,  as  the 
flow  separates  off  the  model  forming  a  shear  layer. 


Figure  72:  Isosurface  of  Q-criteria  around  the  ogive  model  with  starboard  actuation  turned  on  at 
12  m/s. 

To  investigate  the  effect  of  changing  the  strength  of  the  disturbance  created  by  the  forcing 
patches  near  the  tip  of  the  model,  the  moving  wall  velocity  was  varied  from  4  m/s  to  16  m/s  in 
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4  m/s  intervals.  A  high  frequency  component  to  the  side  force  can  be  seen  in  Fig.  71  and  was 
observed  for  the  other  moving  wall  velocities  tested  as  well.  Therefore  to  gauge  the  effect  of 
forcing  on  the  resulting  side  force,  the  final  150  ms  of  each  forcing  cycle  was  used  to  estimate  the 
average  side  force  on  the  model.  Figure  73  shows  the  resulting  time- averaged  side  force  coefficient, 
Cy,  as  a  function  of  the  set  moving  wall  velocity.  To  differentiate  between  port  and  starboard 
forcing,  negative  moving  wall  velocities  indicate  forcing  on  the  port  side,  while  positive  moving 
wall  velocities  indicate  forcing  on  the  starboard  side.  From  Fig.  73,  the  average  side  force  appears 
proportional  to  the  moving  wall  velocity,  at  least  at  this  angle  of  attack  and  Reynolds  number.  This 
same  trend  is  seen  in  the  companion  experiments  when  the  average  side  force  is  compared  to  the 
applied  voltage  to  the  plasma  actuation.Fagley  et  al.  [2012a] 


Figure  73:  Average  side  force  coefficient,  Cy,  for  different  moving  wall  velocities.  A  positive 
moving  wall  velocity  indicates  forcing  on  the  starboard  side,  while  a  negative  moving  wall  veloc¬ 
ity  indicates  forcing  on  the  port  side.  All  forcing  is  tangential  to  the  model  in  the  downstream 
direction. 

5.4.5.4  Proper  Orthogonal  Decomposition  of  the  Flow  Field  The  serial  dataset  alternating 
between  port  and  starboard  forcing  at  12  m/s  (see  Fig.  ??),  which  was  saved  on  the  tap  grid  shown 
in  Fig.  69,  was  analyzed  and  the  spatial/temporal  POD  modes  were  calculated.  Figure  74  shows  the 
cumulative  energy  captured  in  the  calculated  POD  modes  based  on  the  pressure  field.  As  shown  in 
the  inset  of  Fig.  74,  using  the  pressure  field  for  the  POD  analysis  requires  over  50  modes  to  capture 
approximately  99  percent  of  the  energy  in  the  flowfield.  However,  the  very  first  mode  (in  this  case 
the  mean)  captures  98  percent  of  the  energy.  Therefore,  a  large  number  of  modes  are  required  to 
account  for  the  fluctuations  from  the  mean  in  the  pressure  field. 

Figure  75  shows  the  resulting  first  five  spatial  POD  modes  based  on  the  pressure  from  the 
tap  data.  In  this  case,  the  first  mode  exactly  corresponds  to  the  mean  pressure  field,  while  the 
second  spatial  mode  accounts  for  the  shift  of  the  primary  vortices  into  an  asymmetric  state.  As 
shown  in  Fig.  75b),  a  'X'-shaped  pattern  in  the  mode  can  be  seen,  especially  near  the  rear  of  the 
model.  Therefore,  when  the  time  coefficient  is  positive,  the  port  vortex  is  shifted  away  from  the 
ogive,  while  the  starboard  vortex  is  shifted  towards  the  body.  This  shifting  of  the  two  primary 
vortex  positions  is  confirmed  by  the  second  POD  mode  time  coefficient.  Fig.  76,  in  which  the  time 
coefficient  and  the  model  side  force  are  clearly  correlated.  Note  that  in  Fig.  76,  the  time  coefficient 
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Figure  74:  Cumulative  energy  content  captured  in  each  POD  Mode.  The  inset  shows  a  zoomed  in 
view  of  the  first  50  Modes  energy  content. 


has  been  scaled  to  match  the  scaling  of  the  side  force  coefficient.  Furthermore,  very  similar  spatial 
structures  were  also  obtained  when  POD  modes  based  on  total  vorticity  were  calculated. 

As  for  the  higher  spatial  modes  and  associated  time  coefficients,  no  clear  correlation  between 
these  modes  and  the  resultant  side  force  has  been  found.  However,  based  on  the  results  shown 
in  Fig.  74,  these  higher  order  modes  contribute  very  little  to  the  total  amount  of  energy  contained 
in  the  flowfield.  Notice  however  that  these  higher-order  modes  do  not  appear  in  conjugate  pairs 
like  POD  modes  from  periodic  flowfields  which  create  traveling  structures.  This  is  because,  while 
there  are  minor  fluctuations  in  the  positions  of  the  vortex,  once  locked  into  one  state  the  vortices 
tend  to  stay  there. 
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(e) 

Figure  75:  First  five  spatial  POD  modes:  a)  Mode  1,  b)  Mode  2,  c)  Mode  3,  d)  Mode  4,  and  e) 
Mode  5. 
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Figure  76:  Comparison  of  the  instantaneous  side  force  coefficient  and  the  2nd  POD  mode  time 
coefficient. 

As  a  check,  POD  was  performed  on  the  time-resolved  simulation  data.  The  port  and  starboard 
data  sets  at  a  mass  flow  rate  of  0.07  were  concatenated  into  a  single  data  before  performing  the 
analysis.  As  shown  in  Figure  77,  very  similar  results  are  obtained  using  the  time-resolved  data, 
especially  in  the  1st  and  2nd  spatial  modes,  where  the  'X'-shaped  pattern  can  still  clearly  be  seen 
in  the  2nd  mode.  Furthermore,  as  shown  in  Figure  78,  the  2nd  mode  still  closely  follows  the  side 
force  on  the  model.  In  this  case,  a  large  spike  occurs  in  the  3rd  POD  time  coefficient  corresponding 
to  the  swap  in  side  force. 
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(e) 


Figure  77:  First  five  spatial  POD  modes  using  the  time-resolved  simulation  data:  a)  Mode  1,  b) 
Mode  2,  c)  Mode  3,  d)  Mode  4,  and  e)  Mode  5. 


104 


Fagley 


KC  Engineering,  Inc. 


April  8,  2013 


FA7000- 1 0-2-0003 .Final 


Figure  78:  Comparison  of  the  instantaneous  side  force  coefficient  and  the  first  three  POD  mode 
time  coefficient  from  the  time-resolved  data. 
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Begin 

Algorithm 


Figure  79:  Procedure  for  genetic  search  algorithm. 

5.4.6  Sensor  placement 

This  section  presents  a  numerical  study  on  optimal  sensor  location  and  number  of  sensors  to  ade¬ 
quately  estimate  both  the  asymmetric  vortex  states  and  transients  therein.  For  the  purpose  of  this 
study,  the  flow  state  is  defined  as  the  side  force  coefficient  ( Cy )  which  entirely  captures  the  asym¬ 
metric  vortex  behavior.  The  side  force  is  the  desired  control  variable  ,  and  thus  it  is  imperative 
to  have  a  surface  sensor  arrangement  that  accurately  predicts  the  dynamics  of  this  quantity.  The 
following  sections  outlay  the  approach  used  for  determining  surface  sensor  placement  for  the  von 
Karman  ogive. 

5.4.6.1  Optimal  Method  The  optimization  problem  is  solved  by  using  a  constrained,  evolu¬ 
tionary  genetic  algorithm  (GA).  The  constrained  optimization  problem  as  a  stochastic  search  rou¬ 
tine  is  designed  for  the  problem, 

minF(x),  (80) 

xeR« 

where  x  is  the  search  variable  defined  in  space  £  E"  subject  to  an  arbitrary  function  F(x).  This  is 
a  constrained  optimization  problem  with  x  subject  to  the  constraint, 

Li<Xi<Ui ,  (81) 

where  L,-  is  the  lower  bound  and  U\  is  the  upper  bound  for  our  search  vector  xt.  This  evolutionary 
search  based  method  is  a  classical  approach  for  these  multidimensional  constrained  optimizations. 
The  basic  architecture  of  the  GA  is  shown  in  Fig.  79.  With  any  GA  search  routine  an  initial  popu¬ 
lation  or  array  of  potential  solutions  is  randomly  selected.  The  objective  function  is  then  evaluated 
at  each  member  of  the  population.  The  objective  function  provides  a  measure  on  how  well  that 
member  of  the  population  performed.  The  fitness  function  transforms  the  result  into  a  relative 
fitness.  Poorly  fit  members  are  discarded  and  fit  members  are  replicated  in  the  reproduction  step. 
Combinations  of  fit  members  are  randomly  selected  and  then  mutated  to  form  the  next  iterative 
population.  This  process  continues  until  a  desired  stop  criteria  is  satisfied. 
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In  application  to  the  problem  of  the  von  Karman  ogive,  the  objective  function  is  defined  as 

/(xs)  =  Cy(t)  -Cp(xs,t)  (82) 

where  C  is  the  observation  and  computed  from  a  least  squares  estimate  (see  below).  The  fitness 
function  is  defined  as  the  2-norm  of  the  objective  function, 

F(xs)  =  ll/(xs)l|2-  (83) 

The  search  domain  (L,-  <  Xj  <Ui ,  )  is  restrained  to  the  upper  surface  of  the  ogive.  Symmetric 
sensor  placement  is  enforced,  such  that  at  each  axial  position  two  sensors  are  placed  at  the  ±0 
positions. 

5.4.6.2  Linear  Stochastic  Estimation  Linear  stochastic  estimation  (LSE)  is  chosen  as  the 
baseline  estimation  method.  Since  flow  fields  of  interest  are  typically  highly  non-linear,  the  perfor¬ 
mance  of  this  method  usually  tends  to  be  poor  for  most  flow  fields.  However,  a  linear  analysis  is 
always  important  because  it  serves  as  a  benchmark  comparison  for  the  more  complex,  non-linear 
system  ID  methods.  Secondly,  the  computational  time  for  the  linear  analysis  is  negligible,  thus 
allowing  for  an  optimal  study  in  terms  of  sensor  location  and  number. 

The  CFD  simulation  provides  surface  pressure  at  any  point  on  the  surface  of  the  ogive.  The 
array  of  surface  sensors  is  defined  as  the  vector,  xs.  A  linear  mapping  which  estimates  the  state  of 
our  flow  field  is  sought.  The  observation  matrix  C  is  computed  to  best  represent  the  flow  state  in  a 
least  squares  sense.  C  is  computed  by  a  matrix  inversion, 

C  =  Cy(t)p(xs\t)~1  where  C  e  Mmx^.  (84) 

Future  estimates  can  then  be  approximated  by  the  matrix  multiplication, 

Cy(t)  =  Cp(xs,t)label(e.Pest)  (85) 

The  performance  of  the  linear  estimation  approach  is  quantified  by  the  error  norm  of  the  estimate. 
The  error  is  defined  as 

£=  ||Cy(t)-Cp(xs,t)||2  (86) 

This  is  the  very  basic  estimation  method.  The  linear  and  stationary  approach  provides  for  a  means 
to  optimally  solve  for  ideal  surface  arrangement.  A  surface  sensor  array  is  sought  to  minimize 
the  error  in  (86).  Once  the  sensor  array  is  determined,  the  estimation  method  can  be  extended  to 
non- stationary,  higher  order,  non-linear  methods.  The  optimal  arrangement  of  sensors  in  a  linear 
fashion  will  also  be  optimal  for  higher  order  methods  as  well. 

5.4.6.3  Results  The  unforced  data,  shown  in  Fig.  80  a),  indicates  that  while  a  significant  side 
force  coefficient,  \Cy  \  ~  1,  can  be  achieved  with  this  geometry  in  simulations,  it  is  smaller  then 
the  experimental  findings.?  However,  the  unsteadiness  of  the  side  force  coefficient  around  the  zero 
point  shows  that  the  flow  is  unstable  with  respect  to  the  vortex  state  when  no  geometric  disturbance 
is  present  to  enforce  either  the  port  or  starboard  asymmetric  vortex  states. 

Fluidic  actuation  is  introduced  into  the  numeric  simulation  via  two  moving  wall  boundary 
conditions  at  the  tip  of  the  von  Karman  ogive,  as  shown  in  Fig.  68.  The  forcing  is  duty  cycled 
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(a)  Unforced  (b)  Forced 

Figure  80:  a)  Unforced  side  force  coefficient,  b)  Forced  side  force  coefficient. 


between  port  and  starboard  actuators  to  produce  both  asymmetric  states,  deterministically  within 
the  flow  field.  For  both  forcing  situations,  the  unforced  state  of  the  vortices  is  altered  resulting 
in  a  side  force  comparable  to  experimentally  measured  values,  such  that  the  flow  has  entered  into 
a  pseudo-steady  asymmetric  state,  see  Fig.  80  b).  It  was  determined  from  this  data  that  the  data 
for  the  unforced  case  would  provide  more  difficulty  for  the  estimation/optimal  routine  because 
of  the  large  number  of  transients  between  asymmetric  states  as  well  as  the  reduced  magnitude  of 
side  force  (i.e.  smaller  differential  pressures  from  port  to  starboard).  The  unforced  data  was  then 
selected  to  be  the  training  data  for  the  optimal  routine  and  the  forced  data  was  selected  for  the 
validation  of  the  resulting  sensor  configuration. 

The  sensor  location  are  bounded  to  the  leeward  side,  0  =  ±110°,  of  the  ogive  over  an  axial 
range  of  1  <  x/D  <  4.  Another  spatial  constraint  is  the  sensor  array  must  be  symmetric  about  the 
z  =  0  plane  to  ensure  symmetric  sensor  placement.  Thus  the  necessary  parameters  to  optimize 
over  are  the  axial  and  the  azimuthal  positions  of  a  sensor  pair,xs  =  [vi ,  0\,X2,  Oj, . . .  0//T ,  for 

2k  sensors.  The  bounds  or  constraints  on  the  search  criteria  are  thus, 


1  <  x/D  <  4 

o<e<  no0 


(87) 


The  (x,y,z)  coordinates  are  then  computed  from  the  axial  and  theta  orientations,  such  that, 


Xl 

Xl 

VI 

r(x\)cos{6\ ) 

Zl 

r(x\)sin(Q\) 

X2 

x\ 

y2 

r{x\)cos{6\ ) 

.  Z2  . 

r{x\)sin{— (fi) 

(88) 


where  the  radius  at  that  axial  location,  r(x\)  is  computed  from  (??)-  (??).  The  vector  xsurf  contains 
the  coordinate  pairs  over  the  bounded  region  and  mirrored  about  the  meridian  plane.  The  pressure 
information  at  xs  is  linearly  interpolated  from  a  surface  tap  grid  on  the  ogive  using  Delaunay 
triangulation. 
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To  ensure  the  correct  optimal  method  with  associated  parameters  (i.e.  population  size,  mutation 
rate  selection/deletion  properties  and  mutations/crossover  operators)  were  chosen  for  the  optimal 
study  of  sensor  location  and  number  of  sensors,  the  entire  error  surface  was  computed  for  a  single 
sensor  pair  and  compared  to  the  results  from  the  genetic  algorithm.  The  error  surface  plot  is  shown 
in  Fig.  81. 


Figure  81:  Error  surface  of  a  single  sensor  pair  and  optimal  solution  shown  to  be  at  location  of 
minimum  error  thus  validating  the  optimal  routine. 

The  optimal  solution  as  computed  by  the  above  method  produces  xopt/D  =  2.4  and  Qopt  =  1.5 
which  is  plotted  as  a  green  triangle  in  Fig.  81.  This  does  align  with  the  minimum  of  the  error 
surface,  so  the  optimal  method  is  validated  for  this  type  of  minimization  problem.  Note  also  that 
the  error  of  the  linear  estimator  becomes  much  worse  outside  the  bounds  of  x/D  >  3.5  and  x/D  < 
1.5.  This  is  mainly  due  to  the  dynamics  of  the  asymmetric  vortex  phenomena.  The  strongest 
asymmetries  are  seen  in  the  axial  region  from  1.5  >  x/D  <  3.5.  Also  it  is  seen  that  the  minimal 
amount  of  error  is  at  the  largest  azimuthal  positions,  i.e.  6  =  ±110°.  This  indicates  that  the  sensors 
near  the  separation  lines  are  critical  to  capture  the  asymmetric  dynamics. 

The  optimal  solution  was  repeated  for  2, . . . ,  6  pairs  of  sensors.  The  results  are  shown  in  Fig.  82. 
As  the  number  of  sensors  is  increased  the  error  of  the  linear  estimator  is  decreased  for  the  training 
data.  Interestingly  though,  the  error  does  not  decrease  for  the  validation  data  but  rather  increases 
if  more  than  6  sensors  (3  pairs)  are  used.  Because  experimental  implementation  is  crucial  for  this 
project,  it  was  determined  that  a  sensor  array  with  two  sensor  pairs  provided  the  best  compromise. 
Also,  the  prediction  error  was  not  significantly  reduced  when  more  sensor  pairs  were  used.  The 
optimum  positions  for  the  two  sensor  pair  array  were  found  to  be, 

x/D  =[2,3.25] 

e  =  [iios'^s0] 

Figure  83  shows  the  unwrapped  surface  of  the  ogive  colored  in  gray  by  the  mean  pressure  distribu¬ 
tion  for  the  training  data  ensemble.  The  locations  of  the  two  sensor  pairs  are  shown  in  Fig.  83  as 
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Figure  82:  Overall  error  plotted  versus  number  of  sensors  at  their  optimal  locations.  Performance 
of  the  training  and  validation  are  shown  in  blue  and  red,  respectively. 


green  dots.  Figure  83  also  shows  instantaneous  separation  locations  for  the  port  asymmetric  state 
(red)  and  instantaneous  separation  locations  for  the  starboard  asymmetric  state  (blue).  This  figure 
shows  that  the  critical  sensing  locations  for  estimating  the  asymmetric  state  are  near  the  separation 
locations.  Also,  the  separation  locations  do  slightly  move  depending  on  which  asymmetric  state 
is  present.  This  small  fluctuation  in  the  separation  location  provides  for  very  large  fluctuations  in 
relative  pressures  at  the  sensor  location  as  shown  in  Fig.  84  for  each  asymmetric  state. 

The  optimal  routine  found  the  point  at  which  the  largest  differential  pressure  occurred  just  near 
the  separation  location.  The  time  history  of  the  estimation  results  for  the  training  and  validation 
cases  are  shown  in  Fig.  85.  As  shown  the  two  sensor  pairs  predict  the  asymmetric  behavior  for  the 
forced  and  unforced  case  producing  less  than  17%  error  for  each  case.  Some  high  frequencies  are 
not  captured  with  this  sensor  placement,  probably  due  to  the  fact  that  those  asymmetric  pressure 
changes  are  occurring  further  upstream  or  downstream  on  the  body  of  the  ogive. 

5.4.6.4  Experimental  Validation  The  flow  state  estimation  technique  laid  out  in  Fagley  et  al. 
[2012b],  is  experimentally  verified  by  the  following  technique.  The  estimated  side  force,  Cy,  as 
described  by  Eqn.  ??,  needs  to  be  validated  and  compared  to  the  actual  force  on  the  model.  To 
compare  these  two  signals,  the  force  balance  sensor  dynamics  need  to  be  measured  and  modelled. 
For  this  an  impulse  response  to  the  wind  tunnel  model  and  resulting  forces  are  measured  to  model 
the  frequency  response.  Figure  86a  shows  the  frequency  response  of  the  side  force  measurement 
due  to  an  impulse.  A  multi-modal  resonance  is  seen  due  to  the  complex  orientation  of  strain 
gauge/flexure  arrangement  of  the  6  degree  of  freedom  force  balance;  additionally,  each  balance 
channel  shows  a  cross  coupled  behavior  which  is  also  a  factor  for  the  multi-modal  resonance. 
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Figure  83:  Unwrapped  surface  of  the  von  Karman  ogive  colored  by  mean  pressure  distribution. 
Red  lines  indicate  separation  locations  for  port  forcing  and  blue  lines  indicate  separation  location 
for  starboard  forcing.  Optimal  sensor  placement  is  shown  as  green  dots. 


(a)  x/D  =  2  (b)  x/D  =  3.25 

Figure  84:  Sensor  placement  as  well  as  azimuthal  pressure  coefficient  for  port  forcing  (blue)  and 
starboard  forcing  (green)  at  axial  location  a)  x/D  =  2  and  b)  x/D  =  3.25 
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(a)  Training  Data 


(b)  Validation  Data 


Figure  85:  Time  histories  and  estimation  performance  for  training  and  validation  data. 


(b) 

112  Fagley 

Figure  86:  (a)  Frequency  spectrum  of  the  side  force  measurement  due  to  an  impulse  response 
[blue]  and  8th  order  AR  model  [green],  (b)  Time  domain  measurement  of  side  force  coefficient 
[green],  estimated  side  force  from  pressure  signals  [blue]  and  estimated  side  force  measurement 
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Nonetheless  the  spectrum  shown  in  Figure  86b  allows  for  modeling  of  the  sensor  dynamics. 
The  impulse  response  of  this  measurement  device  is  fit  by  an  auto-regressive  system  in  the  form, 

GFB(q)CyFB(t)  =  e(t)  (90) 

where  CyFB  is  the  measurement  of  the  side  force  and  Gfb(<1 )  is  a  frequency  domain  model  of  the 
force  balance.  All  three  signals  can  be  compared  by  the  following  equation, 

GyFB  ~  Gfb(s)  *C  x  P(xs,f).  (91) 

An  example  of  these  signals  due  to  alternating  port  and  starboard  forcing  is  shown  in  Fig¬ 
ure  86b.  The  green  curve  represents  the  actual  side  force  measurement.  The  experimental  mea¬ 
surements  of  the  force  balance  are  compared  directly  to  the  pressure  based  estimates  of  the  force 
as  shown  in  Eq  ??.  Figure  86b  shows  the  actual  force  measurements  in  green,  the  estimated  side 
force  from  the  linear  combination  of  pressure  signals  in  blue,  and  the  estimated  force  measurement 
with  the  force  balance  dynamics  included  in  red.  As  shown  qualitative  agreement  exists  between 
all  three  signals;  the  determination  is  thus  that  the  estimated  side  force  Cy  is  a  more  suitable  signal 
for  the  actual  force  on  the  model,  because  the  sensor  dynamics  of  the  force  balance  are  excluded. 

5.4.6.5  Summary  Unforced  and  forced  CFD  simulations  were  used  to  understand  the  asym¬ 
metric  vortex  state  behavior  on  a  von  Karman  ogive.  The  state  of  the  flow  was  modeled  using  the 
side  force  coefficient.  A  genetic  algorithm  was  used  to  solve  for  optimal  sensor  placement  and 
investigate  the  performance  as  a  function  of  number  of  sensors.  The  fitness  function  of  the  genetic 
algorithm  was  defined  as  the  error  between  the  least  squares  approximation  of  the  surface  pressure 
to  the  defined  flow  state. 

An  optimal  arrangement  of  sensors  was  chosen  which  is  experimentally  feasible  and  in  fact  has 
been  experimentally  implemented.  The  current  study  showed  that  a  total  of  two  sensor  pairs  placed 
at  position  of  x/D  =  [2,3.25]  and  0  =  [±105°,  ±98°],  accurately  predicted  the  flow  state  to  within 
17%  error  for  training  (unforced)  and  validation  (forced)  simulations.  The  placement  algorithm 
showed  that  sensors  placed  very  near  the  separation  point  were  optimal.  This  region  showed  the 
largest  differential  pressures  for  port  and  starboard  asymmetric  vortices,  which  heuristically  is  the 
best  location  for  flow  state  sensors. 
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5.4.7  System  dynamics  and  modeling 

To  model  the  system  dynamics  of  the  asymmetric  vortex  state  behind  the  von  Karman  ogive  at  high 
angles  of  attack,  open-loop  experimental  tests  were  conducted  to  understand  system  characteristics 
in  terms  of  stability/bi- stability,  controllability,  observability,  and  linear/non-linear  behavior.  The 
flow  behind  an  axisymmetric  slender  body  has  previously  been  shown  to  be  completely  bistable  at 
a  sufficiently  large  angle  of  attack  and  small  Reynolds  numbers.  Because  few  tests  have  been  con¬ 
ducted  at  a  Reynolds  number  in  the  range  of  the  current  experiment  (Re  =  156,000),  determining 
if  a  bi-stable  or  proportional  flow  regime  exists  is  critical  in  designing  a  suitable  model  structure 
and  control  system  design. 

The  system  with  the  characteristics  shown  in  Figure  56  is  suitable  for  being  modeled  by  a 
linear  system  with  saturation  points  as  well  as  a  dead  zone.  Therefore,  standard  linear  system 
identification  methods  can  be  used  for  the  system  identification  to  extract  critical  features  such  as 
delay  time,  rise  time,  cut-off  frequency,  phase/gain  margin  and  minimum  phase  behavior.  Once 
these  critical  features  are  well  quantified,  the  linear  system  will  provide  the  means  for  closed- 
loop  controller  development.  For  the  open-loop  database,  the  plasma  actuator  voltage  is  varied 
in  different  manners  to  fully  describe  and  model  the  dynamics.  In  this  experimental  investigation 
three  separate  campaigns  were  conducted  to  develop  the  open-loop  database.  For  training  the 
model,  a  step  response  is  measured,  which  contains  all  necessary  information  for  a  linear  model 
to  be  developed.  Linear  modeling  methods,  such  as  the  output  error,  prediction  error  and  subspace 
identification  methods  are  implemented  to  capture  the  dynamic  response  to  the  step  input.  For 
validation  of  the  developed  model,  both  impulse  and  sinusoidal  forcing  responses  are  compared  to 
experimental  validation. 

Initially,  due  to  geometric  asymmetries,  angular  misalignments  or  flow  imperfections  the  sys¬ 
tem  is  in  the  port  asymmetric  state  which  causes  a  port  attached  vortex,  i.e.  a  negative  Cy  or  ACp. 
This  steady  state  value  is  removed  and  relative  changes  to  the  asymmetric  state  are  analyzed.  For 
all  modeling  purposes  the  side  force  estimate,  Cy,  will  be  used.  Also,  for  all  of  the  data  presented, 
only  the  port  (negative  voltages)  and  starboard  (positive  voltages)  actuators  are  employed  to  influ¬ 
ence  the  flow  state. 

5.4.7.1  Modeling  Techniques  The  system  response  is  modeled  using  a  linear  system  parametriza- 
tion.  The  input  output  relationship  for  this  system  is 

Y(s)  =  Gs(s)U(s).  (92) 

The  structure  of  the  model  in  continuous  time  will  take  the  form, 

r  /  X  N(s)  e^OT  +  qw_i5OT~1  +  aOT-2.yOT~2  •  •  •  a\s  +  ap 
s  D(s)  sn  +  bn-\sn~l  +bn-2Sn~2  ■  ■  -bis  +  bo 

for  a  linear  system  with  m  zeros  and  n  poles  and  a  pure  time  delay,  eGs.  Different  system  identifi¬ 
cation  techniques  exist  for  parameterizing  suitable  orders  G(s)  and  solving  for  coefficients  of  the 
polynomials  in  numerator  and  denominator.  The  three  techniques  for  time  domain  identification 
which  are  examined  in  this  effort  are:  Output  Error  (OE),  Prediction  Error  (PEM),  and  Subspace 
Identification  (SSID)  methods. 
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5.4.7.2  Output  Error  Method  The  output  error  method  is  an  autoregressive  exogenous  input 
(arx)  model  structure,  and  the  identification  method  is  well  representative  of  discrete  time  and 
frequency  domain  data  of  the  form, 

y(k)  =  ^rr\u(k-nk)+e(t).  (94) 

f\Q) 

The  OE  method  minimizes  the  cost  function 

lb(*) -y(fc|©)|ll,  (95) 

given  a  parameterized  vector  which  contains  numerator  order,  denominator  order  and  pure  time 
delay,  rip,  rip,  nk,  respectively.  The  optimal  parameter  vector,  6  is  given  by 

0  =  min^  £  II y(k) -me)\\i  (96) 

For  this  study  the  pure,  convective  time  delay  is  estimated  from  the  step  response  measurements 
as  shown  in  Table  7.  Six  values  for  n k  are  chosen  for  these  formulations.  The  true  pure  delay  time 
is  shown  to  be  8  ms  with  a  sample  time  of  Ts  =  0. 1  ms  which  corresponds  to  a  discrete  delay  time 
of  nk  =  80.  Because  convective  time  and  non-minimum  phase  aren’t  decoupled,  the  convective 
time  parameter  was  varied  to  allow  for  the  zeros  to  adjust  accordingly  to  any  non-minimum  phase 
behavior.  Both  the  numerator  orders,  ns ,  and  denominator  orders, rip,  are  chosen  over  a  range  from 
1  to  5.  With  these  three  parameter  ranges,  a  total  of  150  OE  models  were  compared  and  validated. 


5.4.7.3  Prediction  Error  Method  The  prediction  error  method  (PEM)  has  a  model  structure 
given  by  an  autoregressive  moving  average  (arma)  system,  and  is  an  iterative  identification  ap¬ 
proach  for  multi-input  multi-output  time  domain  data  with  a  model  structure  of  the  form 


A(q)y(k) 


m 

F(q) 


u(k  —  nk )  + 


C(q) 

D(q) 


e(t). 


(97) 


This  linear  time  model  incorporates  a  system  disturbance  term  which  is  filtered  by  the  transfer 
function  (a  type  of  moving  average).  The  parametrization  for  this  model  structure  consisted  of  a 
total  of  six  parameters  as  shown  by, 


0  =  [na,nb,nc,nd,nf,nk,] 


(98) 


where  each  numerator  and  denominator  order  is  denoted  by  nt.  Each  order  was  varied  from  1 
through  5  and  the  delay  term  was  set  to  the  convective  time  delay  computed  from  the  step  input 
which  was  approximately  8  ms.  All  models  were  compared  and  validated  against  experimental 
validation  data  in  following  sections. 


5.4.7.4  Subspace  Identification  Method  The  subspace  identification  (SSID)  method  is  widely 
used  for  black  box  modeling  of  linear  dynamical  systems  directly  in  the  state  space  domain,  which 
are  written  as, 

*k+ 1  =  Axk+Buk  +  Kek, 
yk  =  Cxk  +  Duk  +  ek, 
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where  u *  is  the  m-dimensional  input,  is  the  n-dimensional  state,  is  the  1-dimensional  output 
and  K  is  the  Kalman  gain.  The  SSID  method  is  highly  useful  for  MIMO  systems,  because  of 
its  numerical  robustness,  and  its  model  order  optimization  based  on  the  singular  values  of  the 
Hankel  matrix.  On  the  downside,  large  data  sets  are  needed  to  form  the  block  Hankel  matrices, 
known  deterministic  processes  are  difficult  to  implement,  and  a  strong  theoretical  understanding 
of  observability  and  controllability  is  necessary.  The  basic  premise  is  to  form  the  Hankel  matrices 
of  the  input-output  data  set;  the  observability  matrix, 

0=[C  CA  CA2  •  •  ■CAn~1] T  ,  (100) 

and  the  reversed  controllability  matrix, 

c€=  [An~1B  -A2B  AB  B]T ,  (101) 

are  imbedded  within  this  large  input-output  Hankel  data  matrix.  An  appropriate  model  order  can 
be  estimated  by  the  singular  values  of  this  Hankel  matrix.  Once  the  model  order  is  selected,  system 
matrices  as  shown  in  Eq  (99)  can  be  extracted. 

5.4.7.5  Model  Selection  Experiments  The  three  modeling  techniques  described  in  section 
5.4.7  are  applied  to  the  training  data.  The  RMS  error  between  the  predicted  response  and  ac¬ 
tual  response  serves  as  the  figure  of  merit  for  model  selection.  Transient  areas  where  the  input 
contained  high  frequency  changes  were  weighted  more  heavily  in  the  calculation  of  the  prediction 
error.  The  best  model  from  each  reduced  order  modeling  technique  was  chosen  and  model  order 
and  structure  was  compared.  Model  structure  was  consistent  around  a  4th  or  5th  order  model. 
Also,  the  poles  of  the  models  tended  to  migrate  outside  the  unit  circle  if  the  convective  time  delay 
was  inaccurate;  this  non-minimum  phase  behavior  allowed  for  the  discrepancy  in  the  time  delay  of 
the  model. 

The  results  for  the  best  simulated  model  responses  in  comparison  to  experimental  measure¬ 
ments  of  the  step  response  are  shown  in  Figure  87.  Initial  and  ending  transients  show  very  good 
prediction  of  convective  delay  as  well  as  rise/fall  time  constants  as  shown  in  Table  7.  The  dynam¬ 
ics  vary  slightly  differently  when  the  asymmetric  state  transitions  from  port  to  starboard  versus 
from  starboard  to  port.  The  transient  time  from  starboard  to  port  vortex  states  is  shorter  and  the 
overshoot  is  greater,  as  Table  7  indicates.  This  is  potentially  because  the  initial  state  prefers  the 
port  asymmetric  state  which  may  provide  an  additional  restoring  force  to  the  vortex  dynamics. 
The  linear  approach  taken  in  this  paper  finds  the  mean  dynamics  between  each  state  trajectory  as 
shown  by  Figure  87.  Nevertheless,  each  of  the  models  replicates  the  asymmetric  vortex  dynamics 
to  a  step  input  very  well,  thus  validating  the  model  parameterizations  as  well  as  model  selection 
technique. 
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Figure  87:  Step  response  comparison  of  OE,  PEM,  and  SSID  modeling  methods,  (a)  Initial  tran¬ 
sient  and  (b)  ending  transient  of  step  response  of  estimated  side  force,  Cy,  based  on  pressure  mea¬ 
surements 

The  validation  data  sets  were  also  used  to  evaluate  model  performance.  The  models  were 
calculated  against  both  of  the  sinusoidal  forcing  and  impulse  forcing  inputs.  The  response  of 
the  asymmetric  state  and  model  responses  were  compared.  Figure  88a  shows  the  summary  of 
the  frequency  response  data  which  aligns  well  with  the  raw  frequency  measurements.  The  phase 
relationship  is  also  shown  in  Figure  88b.  To  select  between  the  three  different  model  development 
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approaches,  the  error  is  minimized  in  the  frequency  domain.  The  prediction  error  technique  most 
adequately  fits  the  frequency  domain  data,  in  both  magnitude  and  phase. 

Interestingly,  the  cutoff  frequency  of  the  system  which  is  determined  from  a  —90 deg  phase 
is  approximately  1/(2t).  This  means  more  or  less  that  any  frequencies  larger  than  an  associated 
period  of  two  flow  through  times  will  be  greatly  attenuated.  This  is  shown  in  Figure  88a. 


(a) 


(b) 

Figure  88:  Comparison  of  model  frequency  response  and  experimental  measurements. 


Figure  89  shows  the  impulse  response  with  varying  duty  cycles  for  the  experimental  and  sim- 
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ulated  output  error  model  results.  The  various  colors  represent  the  different  duty  cycles.  The 
solid  lines  represent  experimental  measurements  and  the  dashed  lines  represent  the  PEM  model 
prediction.  All  of  the  impulses  were  initiated  at  time  equal  to  zero  with  the  ending  duration  of 
the  impulse  indicated  by  a  vertical  dashed  line.  The  linear  model  has  a  slightly  different  gradient 
during  transient  times  and  a  small  amount  of  overshoot  when  returning  to  the  initial  state. 


Figure  89:  Validation  of  output  error  model  for  simulation  of  impulse  response  with  varying  duty 
cycles.  Model  response  is  shown  in  dashed  line  and  experimental  measurement  is  shown  in  solid 
line 

5.4.7.6  Model  Selection  Simulations  The  three  modeling  techniques  described  in  section  5 .4.7 
are  applied  to  the  step  response  data.  The  RMS  error  between  the  predicted  response  and  actual 
response  serves  as  the  figure  of  merit  for  model  selection.  Transient  areas  where  the  input  contained 
high  frequency  changes  were  weighted  more  heavily  in  the  calculation  of  the  prediction  error.  The 
best  model  from  each  reduced  order  modeling  technique  was  chosen  and  model  order  and  structure 
was  compared.  Model  structure  was  consistent  around  a  4th  or  5th  order  model.  Also,  the  poles 
of  the  models  tended  to  migrate  outside  the  unit  circle  if  the  convective  time  delay  was  inaccurate; 
this  non-minimum  phase  behavior  allowed  for  the  discrepancy  in  the  time  delay  of  the  model. 

The  results  for  the  best  simulated  model  responses  in  comparison  to  experimental  measure¬ 
ments  of  the  step  response  are  shown  in  Figure  90.  Initial  and  ending  transients  show  very  good 
prediction  of  convective  delay  as  well  as  rise/fall  time  constants  as  shown  in  Table  7.  The  dynam¬ 
ics  vary  slightly  differently  when  the  asymmetric  state  transitions  from  port  to  starboard  versus 
from  starboard  to  port.  The  transient  time  from  starboard  to  port  vortex  states  is  shorter  and  the 
overshoot  is  greater,  as  Table  8  indicates.  This  is  potentially  because  the  initial  state  prefers  the 
port  asymmetric  state  which  may  provide  an  additional  restoring  force  to  the  vortex  dynamics. 
The  linear  approach  taken  in  this  paper  finds  the  mean  dynamics  between  each  state  trajectory  as 
shown  by  Figure  90.  Nevertheless,  each  of  the  models  replicates  the  asymmetric  vortex  dynamics 
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to  a  step  input  very  well,  thus  validating  the  model  parameterizations  as  well  as  model  selection 
technique. 


Figure  90:  Step  response  comparison  of  OE,  PEM,  and  SSID  modeling  methods,  (a)  Initial  tran¬ 
sient  and  (b)  ending  transient  of  step  response  of  estimated  side  force,  Cy,  based  on  pressure  mea¬ 
surements 

The  validation  data  sets  were  also  used  to  evaluate  model  performance.  The  models  were 
calculated  against  both  of  the  sinusoidal  forcing  and  impulse  forcing  inputs.  The  response  of 
the  asymmetric  state  and  model  responses  were  compared.  Figure  91a  shows  the  summary  of 
the  frequency  response  data  which  aligns  well  with  the  raw  frequency  measurements.  The  phase 
relationship  is  also  shown  in  Figure  91b.  To  select  between  the  three  different  model  development 
approaches,  the  error  is  minimized  in  the  frequency  domain.  The  prediction  error  technique  most 
adequately  fits  the  frequency  domain  data,  in  both  magnitude  and  phase. 

Interestingly,  the  cutoff  frequency  of  the  system  which  is  determined  from  a  —90  deg  phase 
is  approximately  1  /  (2t).  This  means  more  or  less  that  any  frequencies  larger  than  an  associated 
period  of  two  flow  through  times  will  be  greatly  attenuated.  This  is  shown  in  Figure  91a. 
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Figure  91:  Comparison  of  model  frequency  response,  both  magnitude  (a)  and  phase  (b). 

Figure  92  shows  the  impulse  response  with  varying  duty  cycles  for  the  CFD  simulation  results 
and  simulated  prediction  error  model  results.  The  various  colors  represent  the  different  duty  cy¬ 
cles.  The  solid  lines  represent  experimental  measurements  and  the  dashed  lines  represent  the  PEM 
model  prediction.  All  of  the  impulses  were  initiated  at  time  equal  to  zero  with  the  ending  duration 
of  the  impulse  indicated  by  a  vertical  dashed  line.  The  linear  model  has  a  slightly  different  gradient 
during  transient  times  and  a  small  amount  of  overshoot  when  returning  to  the  initial  state. 
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Figure  92:  Validation  of  prediction  error  model  for  simulation  of  impulse  response  with  varying 
duty  cycles.  Model  response  is  shown  in  dashed  line  and  CFD  simulation  results  are  shown  in 
solid  line 
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5.4.8  Closed  loop  dynamics 

Now  that  the  system  dynamics  in  Eq  92  have  been  modeled,  i.e.  the  relationship  between  plasma 
voltage  and  estimated  side  force  response  is  determined,  the  closed-loop  system  can  be  realized. 
The  overall  design  of  the  control  diagram  is  shown  in  Figure  93.  The  type  of  control  system 
selected  is  a  reference  tracking  feed  forward  approach  where  Gs(s)  is  the  plant  as  developed  in 
section  5.4.7,  Gc(s)  is  the  control  system,  and  G(i(s)  is  the  output  disturbance/measurement  noise. 
The  unforced  fluid  state  and  measurement  noise  is  modeled  by  an  output  disturbance,  which  is 
colored  by  the  unforced  dynamics  of  the  sensor  measurements.  The  output  disturbance,  Gd(s), 
may  be  represented  by  the  unforced,  natural  fluctuating  state  of  the  flow.  An  autoregressive  model 
is  suitable  for  the  determination  of  this  system.  The  success  of  the  feedback  control  scheme  will 
be  determined  if  adequate  disturbance  rejection  as  well  as  reference  tracking  ability  are  shown. 


Figure  93:  Closed-loop  block  diagram  of  output  disturbance  model  for  controller  verification. 

The  closed-loop  system  is  formulated  such  that 

y  GSGC  Gd 

[  1 +GsGd  1 +GsGd 

where  r  and  d  are  the  reference  and  disturbance  inputs,  respectively.  The  transfer  function  between 
different  input-output  pairs  can  be  analyzed  for  varying  forms  of  Gc(s).  For  the  purpose  of  this 
paper,  the  design  of  the  controller,  Gc(s),  is  standard  PID  control.  A  PID  control  algorithm  is 
implemented  because  of  the  simplicity  and  ease  of  design.  The  asymmetric  vortex  dynamics  lend 
themselves  very  well  to  linear  time  invariant  systems,  so  a  simple  control  algorithm  is  appropriate 
for  control  of  the  vortex  flow  behind  the  ogive.  The  control  algorithm  is  given  by 

Gc(s)=Kp  +  Kds  +  1j,  (103) 

where  Kp,Kd,  Ki  are  the  proportional,  derivative  and  integral  terms,  respectively.  A  standard  tuning 
method  is  adopted  where  the  gains  are  varied  in  a  systematic  fashion  to  achieve  proper  closed-loop 
response  to  a  step  reference  input. 

5.4.8.1  Closed  Experimental  Model  Results  The  response  of  GSGC/(1  +  GsGd )  is  shown  for 
varying  proportional  and  integrator  gains  in  Figure  94.  Selected  gains  for  the  PI  controller  are, 
Ki  =  80  and  Kp  =  1 .2.  The  derivative  term  caused  an  instability  in  the  transfer  function,  GSGC/ (1  + 
GsGd ),  purely  due  to  the  time  delay  in  the  system. 


r 

d 


(102) 
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Figure  94:  (a)  Closed-loop  step  response  with  varying  proportional  gain,  (b)  Closed-loop  step 
response  with  varying  integral  gain. 
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Figure  95:  Magnitude  (a)  and  phase  (b)  response  of  open  loop  system,  G0i  =  GSGC.  The  gain  and 
phase  margins  are  computed  as  2.44 dB  and  62°,  respectively. 

The  frequency  response  for  the  reference  tracking  and  disturbance  rejection  capabilities  are 
shown  in  Figure  96.  As  shown  in  Figure  96a,  the  closed-loop  system  response  adequately  follows 
reference  signals  up  to  approximately  1  / (It)  which  was  determined  to  be  the  cutoff  frequency 
in  the  open-loop  analysis  of  the  dynamics.  Figure  96b  shows  the  closed-loop  system  attenuation 
of  disturbances  (i.e.  the  ability  of  the  closed-loop  system  to  reduce  fluctuations  as  a  function  of 
frequency).  Disturbances  are  attenuated  up  to  a  frequency  of  1/ (4t)  which  turns  out  to  be  half  of 


125 


Fagley 


KC  Engineering,  Inc. 


April  8,  2013 


FA7000- 1 0-2-0003 .Final 


the  cut  off  frequency. 


Figure  96:  (a)  Closed-loop  system  response  for  input  reference  to  output  response,  (b)  Closed-loop 
system  response  for  input  disturbance  to  output  response. 

A  typical  time  simulation  is  shown  in  Figure  98  to  a  time  varying  reference  with  a  uniform 
random  disturbance  input.  As  shown  the  response  of  the  side  force  adequately  follows  a  reference 
signal.  The  controller  is  designed  aggressively  enough  to  have  over/under  shoot  characteristics 
with  step  changes.  Additionally,  the  disturbances  at  lower  frequencies  are  reduced  in  size. 
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Figure  97 :  Time  simulation  of  closed-loop  system  for  varying  reference  and  disturbance  excitation. 

5.4.9  Closed  Experimental  Results 

Once  the  closed-loop  model  was  fully  developed  and  analyzed  (as  discussed  in  Section  ??)  a  ref¬ 
erence  tracking  exercise  was  conducted  with  the  experimental  model  in  the  wind  tunnel  to  asses 
the  performance  of  the  closed-loop  controller.  For  these  investigations  an  arbitrary  piecewise  refer¬ 
ence  waveform  was  generated  for  the  controller  to  track,  where  the  target  side  force  coefficient  was 
changed  seven  times  across  the  15  sec.  testing  period  to  both  positive  and  negative  side  force  coef¬ 
ficients.  A  PID  (Proportional-Integral-Derivative)  controller  was  developed  utilizing  both  the  port 
and  starboard  actuators  to  impart  control  and  all  four  time-resolved  pressure  transducers  were  used 
to  estimate  the  instantaneous  side  force  coefficient  on  the  body  for  feedback.  Figure  98  displays 
the  performance  of  the  experiment  for  the  reference  tracking  exercise,  where  the  experimental  sig¬ 
nal  is  a  phase  averaged  result  of  five  independent  experimental  tests  following  the  same  arbitrary 
reference  signal  (also  presented).  For  the  experiments  the  proportional  and  integral  gains  were 
0.25  and  0.000977,  respectively.  The  derivative  gain  was  set  to  zero  because  it  was  found  during 
the  modeling  that  any  amount  of  derivative  gain  forced  the  model  to  go  unstable.  Clearly,  the  ex¬ 
periment  was  successfully  able  to  track  the  reference  signal  in  the  mean  of  the  linearly  estimated 
side  force  coefficient,  however  significant  higher  frequency  fluctuations  were  still  observed,  which 
the  controller  was  unable  to  modify. 
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Figure  98:  Closed  loop  experimentation  results  with  same  controller  from  above  section. 

5.4.9.1  Closed  Simulation  Model  Results  The  response  of  GsGc/(  1  +GsGd)  is  shown  for 
varying  proportional  and  integrator  gains  in  Figure  94.  Selected  gains  for  the  PI  controller  are, 
Ki  =  .012  and  Kp  =  l  x  10  7 .  The  derivative  term  caused  an  instability  in  the  transfer  function, 
GsGc/(  1  +  GSGC(),  purely  due  to  the  time  delay  in  the  system. 


Figure  99:  (a)  Closed-loop  step  response  with  varying  proportional  gain,  (b)  Closed-loop  step 
response  with  varying  integral  gain. 


128 


Fagley 


KC  Engineering,  Inc. 


April  8,  2013 


FA7000- 1 0-2-0003 .Final 


Figure  100:  Magnitude  (a)  and  phase  (b)  response  of  open  loop  system,  G0i  =  GSGC.  The  gain  and 
phase  margins  are  computed  as  U dB  and  67°,  respectively. 

The  frequency  response  for  the  reference  tracking  and  disturbance  rejection  capabilities  are 
shown  in  Figure  101.  As  shown  in  Figure  101a,  the  closed-loop  system  response  adequately  fol¬ 
lows  reference  signals  up  to  approximately  1  /  (2t)  which  was  determined  to  be  the  cutoff  fre¬ 
quency  in  the  open-loop  analysis  of  the  dynamics.  Figure  101b  shows  the  closed-loop  system 
attenuation  of  disturbances  (i.e.  the  ability  of  the  closed-loop  system  to  reduce  fluctuations  as  a 
function  of  frequency).  Disturbances  are  attenuated  up  to  a  frequency  of  1/ (4t)  which  turns  out  to 
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be  half  of  the  cut  off  frequency. 


Figure  101:  (a)  Closed-loop  system  response  for  input  reference  to  output  response,  (b)  Closed- 
loop  system  response  for  input  disturbance  to  output  response. 

A  typical  time  simulation  is  shown  in  Figure  98  to  a  time  varying  reference  with  a  uniform 
random  disturbance  input.  As  shown  the  response  of  the  side  force  adequately  follows  a  reference 
signal.  The  controller  is  designed  aggressively  enough  to  have  over/under  shoot  characteristics 
with  step  changes.  Additionally,  the  disturbances  at  lower  frequencies  are  reduced  in  size. 
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Figure  102:  Time  simulation  of  closed-loop  system  for  varying  reference  and  disturbance  excita¬ 
tion. 

5.4.9.2  Closed  Simulation  Results  Using  the  unforced  data,  a  linear  prediction-error  mini¬ 
mization  method  was  used  to  model  the  dynamics  of  the  flowfield  for  different  momentum  coef¬ 
ficients.  Based  on  the  model,  a  PID  controller  was  developed  to  reference  track  a  prescribed  side 
force  trajectory.  The  details  of  the  model  development  as  well  as  the  PID  controller  can  be  found 
in  Fagley  et  al.  [submitted  2012]  (experimental)  as  well  as  Porter  et  al.  [2013]  (computational). 
Figure  103  shows  the  initial  results  obtained  using  the  PID  controller  in  conjunction  with  Cobalt. 
A  reference  side  force  of  Cy  =  0.5  was  used  to  test  if  the  controller  could  reference  track.  For  this 
simulation,  an  aggressive  gain  for  the  proportional  and  integral  components  was  used.  As  a  result, 
the  side  force  overshoots  its  reference  condition.  At  this  point,  the  controller  turns  on  the  starboard 
actuator  to  counteract  this  overshoot,  creating  a  large  oscillation  in  the  side  force.  However,  note 
that  this  is  exactly  what  was  predicted  in  the  model  simulation  of  the  controller  (Fig.  ??).  While 
there  is  a  small  discrepancy  in  the  transition  time  between  the  controller  model  and  the  CFD  simu¬ 
lation,  the  general  overall  trends  are  captured  in  the  model.  As  shown  in  the  model,  this  oscillation 
from  aggressive  PID  gains  eventually  damps  out  and  the  controller  is  able  to  stabilize  the  side  force 
at  the  desired  reference.  It  is  postulated  that  if  the  current  CFD  simulation  were  carried  out  farther, 
the  same  results  would  be  seen,  especially  since  the  overshoot  of  the  second  peak  in  the  CFD  is 
smaller  then  the  initial  overshoot,  indicating  that  it  is  starting  to  be  damped  out. 

5.4.10  Modeling  summary 

The  asymmetric  vortex  regime  of  a  von  Karman  ogive  with  a  fineness  ratio  of  3.5  is  experimentally 
studied  at  a  Reynolds  number  of  156,000.  Both  port  and  starboard  plasma  actuators  are  used  to 
introduce  fluidic  disturbances  at  the  tip  of  the  ogive.  These  disturbances  are  amplified  through  the 
flow’s  convective  instability  to  produce  a  deterministic  port  or  starboard  asymmetric  vortex  state 
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Figure  103:  Closed-loop  simulation  with  a  reference  Cy  =  0.5  and  the  corresponding  out  of  the 

port  and  starboard  actuators. 

(i.e.  side  force).  Accurate  control  or  manipulation  of  this  asymmetric  vortex  phenomenon  holds 
the  potential  for  increased  maneuverability  and  stability  characteristics  of  slender  flight  vehicles. 

Unforced  and  open-loop  experimental  tests  were  carried  out  to  understand  and  quantify  the  vor¬ 
tex  dynamics.  Step,  impulse  and  sinusoidal  modulation  inputs  provided  the  necessary  dynamics 
and  diverse  training  and  validation  data  sets  for  the  formulation  of  a  linear  time  invariant  dynamical 
model.  Standard  linear  system  identification  approaches  were  implemented  to  represent  the  train¬ 
ing  data  set.  In  particular,  output  error,  prediction  error  and  subspace  identification  methods  were 
used  to  capture  the  asymmetric  vortex  dynamics.  These  methods  were  validated  by  time  and  fre¬ 
quency  domain  methods.  The  measurements  and  modeling  methods  showed  the  cutoff  frequency 
of  the  flow  to  be  around  50  Hz  which  is  directly  related  to  two  flow  through  times,  i.e.  the  time  it 
takes  a  particle  to  flow  from  the  tip  of  the  model  to  the  base  of  the  ogive  section. 

A  closed-loop  system  was  designed  such  that  the  unforced  fluid  dynamics  and  measurement 
noise  were  modeled  as  an  output  disturbance.  The  prediction  error  model  was  well  suited  for  this 
system.  A  PID  controller  was  implemented  in  the  closed  loop  system  and  designed  for  adequate 
disturbance  rejection  and  reference  tracking  performance.  The  closed  loop  transfer  functions  were 
analyzed.  A  time  simulation  was  shown  in  which  the  controller  was  able  to  guide  the  asymmet¬ 
ric  vortex  state  to  an  arbitrary  asymmetric  pressure  distribution  while  adequately  regulating  the 
disturbances.  To  improve  this  control  design  approach  a  predictor  model  would  be  essential  to 
reduce  the  convective  time  delay  from  the  actuator  to  the  sensor.  Alternatively,  the  sensors  would 
have  to  be  placed  closer  to  the  nose  of  the  ogive  which  would  reduce  the  amplitude  of  the  pressure 
measurements,  reducing  the  signal  to  noise  ratio. 
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5.5  Business  Summary 

The  business  plan  was  divided  into  both  a  salary  portion  and  travel  portion  for  the  senior  engineer. 
The  projected  salary  portion  is  listed  in  table  9,  itemized  yearly  for  the  duration  of  the  contractual 
agreement.  The  projected  salary  portion  is  listed  in  10  for  travel  to  one  conference  per  year. 


Table  9:  Total  projected  salary  for  contract  duration 


Period  of  Performance 

Labor  Category 

Rate 

Hours 

Yearly  Price 

8  Jan  2010  to  7  Jan  201 1 

Senior  Engineer 

$35.51 

1880 

$66,750 

8  Jan  201 1  to  7  Jan  2012 

Senior  Engineer 

$39.89 

1880 

$75,000 

8  Jan  2012  to  7  Jan  2013 

Senior  Engineer 

$42.55 

1880 

$80,000 

Total 

$221,750 

Table  10:  Total  projected  travel  for  contract  duration 


Travel  Expenses 

Flights 

Lodging/night 

Per  Diem  (M&IE) 

Registration 

Total 

2010  Conference 

$800.00 

$120.00 

$80.00 

$300.00 

$2,100.00 

2011  Conference 

$800.00 

$120.00 

$80.00 

$300.00 

$2,100.00 

2012  Conference 
Total 

$800.00 

$120.00 

$80.00 

$300.00 

$2,100.00 

$6,300.00 

The  actual  costs  incurred  are  listed  in  Tables  11  and  12  for  salary  and  travel,  respectively.  As 
shown,  costs  between  projected  and  actual  did  differ  slightly.  This  was  mainly  due  to  obligatory 
issues  during  Y3.Q1  and  Y3.Q2.  Also,  more  travel  was  required  over  the  course  of  the  contractual 
duration. 
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Table  11:  Total  salary  paid  from  January  2010  to  January  2013 


Period  of  Performance  Rate 

Hours 

Price 

Y1.Q1 

$35.51 

470 

$  16,689.70 

Y1.Q2 

$35.51 

470 

$  16,689.70 

Y1.Q3 

$35.51 

470 

$  16,689.70 

Y1.Q4 

$35.51 

470 

$  16,689.70 

Y2.Q1 

$  39.90 

470 

$  18,753.00 

Y2.Q2 

$  39.90 

470 

$  18,753.00 

Y2.Q3 

$  39.90 

470 

$  18,753.00 

Y2.Q4 

$  39.90 

470 

$  18,753.00 

Y3.Q1 

$  43.67 

450 

$  19,651.50 

Y3.Q2 

$  43.67 

420 

$  18,341.40 

Y3.Q3 

$  43.67 

470 

$  20,524.90 

Y3.Q4 

$  43.67 

470 

$  20,524.90 

Total  Salary  Paid 

$  220,813.50 

Table  12:  Total  salary  paid  from  January  2010  to  January  2013 


Dates 

Type 

Location 

Total 

May  26-28,  2010 

AFC  II 

Germany 

$  1,575.00 

August  1-5,  2010 

AIAA  GNC 

Toronto 

$  2,053.76 

Dec  7-10,  2011 

ADD  CoOp 

South  Korea 

$2,100.00 

Dec  2-7,  2012 

ADD  CoOp 

South  Korea 

$2,100.00 

Total  Travel  Funds 

$  8,028.76 
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